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:
Hi Erwin,
ReplyDeleteCould you explain how is made to install the subsystem "nls" in GAMS?
Thanks and congratulations for your great blog, it is very useful!
See http://www.amsterdamoptimization.com/statistics.html.
ReplyDeleteHi, Erwin, great code!
ReplyDeletecan 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!!!
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.
ReplyDeletethanks, Erwin, you are right. sorry it should be
ReplyDeletesolve 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!
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?
ReplyDeleteHi Utkur,
DeleteDid you get any answer for your question? If yes, I would be interested ...
Thanks
Marie