Monday, November 23, 2009

Stochastic Forecast with GAMS

A client of mine asked me to see if we could implement some stochastic forecasting scheme in GAMS and incorporate this into their model. A SAS program to calculate stochastic forecasts was provided, and it was easy to translate this into GAMS:

sets
   t    
/1993*2050/
   tr(t)
'used in regression' /1993*2009/
   tf(t)
'used in forecast' /2010*2050/
;

parameter yield(t) /

  
1993    4.166666667
  
1994    4.454901961
  
1995    4.111111111
  
1996    4.558823529
  
1997    6.097637795
  
1998    5.182341651
  
1999    5.548387097
  
2000    5.464868701
  
2001    6
  
2002    6.326530612
  
2003    6.52173913
  
2004    7.374100719
  
2005    6.475409836
  
2006    8.035714286
  
2007    6.395705511
  
2008    7.003890813
  
2009    7.339368303

/;

parameter year(t);
year(t) = 1992+
ord(t);

*-----------------------------------------------------------------------------
* regression
*-----------------------------------------------------------------------------

variables
   sse
   b0
'constant term'
   b1
;
equations
   dummy
   fit(tr)
;

dummy.. sse =n= 0;
fit(tr).. yield(tr) =e= b0 + b1*year(tr);

option lp=ls;
model ols /dummy,fit/;
solve ols using lp minimizing sse;

* get standard error of regression
* aka Root MSE in SAS
scalars rmse;
execute_load 'ls.gdx',rmse=sigma;
display rmse;

*-----------------------------------------------------------------------------
* forecast
*-----------------------------------------------------------------------------

* this if you want every run to be different
* execseed = 1+gmillisec(jnow);

parameter
   statfcst(t)  
'static forecast'
   stochfcst(t) 
'stochastic forecast'
;
statfcst(tr) = yield(tr);
statfcst(tf) = b0.l + b1.l*year(tf);

stochfcst(tr) = yield(tr);
stochfcst(tf) = b0.l + b1.l*year(tf) + rmse*normal(0,1);

display statfcst,stochfcst;

*-----------------------------------------------------------------------------
* Monte Carlo Simulation
*-----------------------------------------------------------------------------

set r 'replications' /1*10000/;
parameter allyield(r,tf);
allyield(r,tf) = b0.l + b1.l*year(tf) + rmse*normal(0,1);

parameter
   stochfcst2(t) 
'stochastic forecast'
;

stochfcst2(tr) = yield(tr);
stochfcst2(tf) =
sum(r, allyield(r,tf))/card(r);
display stochfcst2;

It is not very difficult to incorporate a SOLVE in the last part. We need a real loop over r for that. The code can look like:

parameter result(r,t) 'used to collect model solutions';
loop(r,
   yield(tf) = b0.l + b1.l*year(tf) + rmse*normal(0,1);
* do something with yield(t) in model
  
solve gams_model using nlp minimizing z;
* collect solutions and yields
   allyield(r,t) = yield(t);
   result(r,t) = x.l(t);
);
* get means
parameter myield(t),mresult(t);
myield(t) =
sum(r,allyield(r,t))/card(r);
mresult(t) =
sum(r,result(r,t))/card(r);


GAMS has some cool looking fan-charts, may be I can use this in this project!




References:

7 comments:

  1. Hi Erwin,

    Could you explain how is made to install the subsystem "nls" in GAMS?

    Thanks and congratulations for your great blog, it is very useful!

    ReplyDelete
  2. Hi, Erwin, great code!

    can i ask what if r (the replication) is an index of a variable rather than a parameter?

    for example, say allyield (r,tf) is a variable of the equation OPT, and we are going to find the optimal solution of allyield(r,tf) for every r.

    i tried to write something, but i guess it is wrong:

    SET r 'replications' /1*10000/;
    VARIABLE allyield(r,tf)
    z(r);
    EQUATION define(r,tf)
    OPT(r);
    define(r,tf).. allyield(r,tf) =E= b0.l + b1.l*year(tf) + rmse*normal(0,1)+ a(tf);
    * a(tf) is another variable
    OPT(r).. z(r) =E= allyield (r,tf) + b(tf);
    * b(tf) is another variable

    MODEL MDL /define, opt/;
    LOOP (r,
    solve MDL using lp minimizing OPT(r)
    );

    Do you think the code like above can work? I think it's very unlikely because I've never read any code written with an index attached after the object variable (here is OPT)...it seems to be an incorrect statement.

    What's your opinion? How would you correct it?

    GREAT THANKS!!!

    ReplyDelete
  3. No that does not work. There are at least a few things wrong with your approach (OPT is an equation, and the objective variable should be a scalar variable). I suggest to contact support@gams.com to get help on your particular model.

    ReplyDelete
  4. thanks, Erwin, you are right. sorry it should be

    solve MDL using lp minimizing z(r);

    but again this ain't gonna work 'cause theres still an index attached.. but well i've got a solution =) just to throw off all the index r from all the variables, and only include it in the loop body.

    Thanks anyway!

    ReplyDelete
  5. Hi Erwin, I have a question. How to include several uncertainty aspects into the model? For instance, in my farm model I want to do monte carlo simulation for irrigation water, price, and yield variability. Could you help how in this case loop command will look like?

    ReplyDelete
    Replies
    1. Hi Utkur,

      Did you get any answer for your question? If yes, I would be interested ...

      Thanks
      Marie

      Delete