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/

;

parameteryield(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

/;

parameteryear(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);

optionlp=ls;modelols /dummy,fit/;solveols using lp minimizing sse;

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

*-----------------------------------------------------------------------------* 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);

displaystatfcst,stochfcst;

*-----------------------------------------------------------------------------* Monte Carlo Simulation*-----------------------------------------------------------------------------setr 'replications' /1*10000/;parameterallyield(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);displaystochfcst2;

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:

parameterresult(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

solvegams_model using nlp minimizing z;* collect solutions and yields

allyield(r,t) = yield(t);

result(r,t) = x.l(t);

);* get meansparametermyield(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