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