## 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;
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:

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!

2. 3. 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!!!

4. 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.

5. 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!

6. 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?

1. Hi Utkur,

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

Thanks
Marie