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: