Conditional Forecasts and Scenario Analysis
/*--------------------------------------------------------------
SAS Sample Library
Name: varex05.sas
Description: Example program from SAS/ETS User's Guide,
The VARMAX Procedure
Title: Conditional Forecasts and Scenario Analysis
Product: SAS/ETS Software
Keys: Vector AutoRegressive Moving-Average
processes with eXogenous regressors
PROC: VARMAX
Notes:
--------------------------------------------------------------*/
****************************************************************;
* Hard Conditions ;
****************************************************************;
title 'Conditional Forecasts and Scenario Analysis';
%macro cfSimulateData(dgpi,T,lead,tblDGP,tblSample);
* dgpi: index of DGP;
* T: in-sample sample size;
* lead: future horizons;
* tblDGP: output table name for full-sample data;
* tblSample: output table name for in-sample data;
proc iml;
* simulate the data;
* trivariate VAR(2) model;
seed = 12345 + &dgpi.; * random seed;
n = 3; * dim of dependent variable;
T = &T.; * in-sample sample size;
lead = &lead.; * future horizons;
p = 2; * AR order;
* parameter values;
const = {2.425, 0.054, -0.110};
phi = {0.234 -0.134 -0.057,
0.029 0.575 0.200,
0.059 0.038 1.006,
0.164 -0.150 -0.165,
-0.039 0.138 -0.184,
0.031 0.019 -0.087};
sigma = {9.265 0.296 0.553,
0.296 1.746 0.184,
0.553 0.184 0.752};
call varmasim(y,phi) sigma = sigma n = T+lead seed = seed;
mu = (inv(I(3)-phi[1:3,]-phi[4:6,])*const)`;
y = y + mu;
name={y1 y2 y3};
create &tblDGP. from y[colname=name];
append from y;
close;
quit;
data &tblSample.; set &tblDGP.(obs=&T.); run;
%mend;
%macro hcConstructScenarios(T,tblDGP,tblScenarios);
* T: in-sample sample size;
* tblDGP: input table name for full-sample data;
* tblScenarios: output table name for scenarios;
data &tblScenarios.;
set &tblDGP.(firstobs=%eval(&T.+1) obs=%eval(&T.+1))
&tblDGP.(firstobs=%eval(&T.+1) obs=%eval(&T.+2))
&tblDGP.(firstobs=%eval(&T.+1) obs=%eval(&T.+3))
&tblDGP.(firstobs=%eval(&T.+1) obs=%eval(&T.+4))
&tblDGP.(firstobs=%eval(&T.+1) obs=%eval(&T.+1));
* scenario 1: y3(1) is known;
if(_N_=1) then do;
y1=.; y2=.; myscenario=1;
end;
* scenario 2: y3(1:2) is known;
if(_N_>1 and _N_<=3) then do;
y1=.; y2=.; myscenario=2;
end;
* scenario 3: y3(1:3) is known;
if(_N_>3 and _N_<=6) then do;
y1=.; y2=.; myscenario=3;
end;
* scenario 4: y3(1:4) is known;
if(_N_>6 and _N_<=10) then do;
y1=.; y2=.; myscenario=4;
end;
* scenario 5: nothing is known (unconditional forecast);
if(_N_>10) then do;
y1=.; y2=.; y3=.; myscenario=5;
end;
run;
%mend;
%macro hcEstimateAndForecast(tblSample,tblScenarios,alpha,lead,nmc,
tblF,tblCf,tblCfSim);
* tblSample: input table name for in-sample data;
* tblScenarios: input table name for scenarios;
* alpha: size of the confidence interval or credible interval;
* lead: future horizons;
* nmc: number of Monte Carlo iterations;
* tblF: output table name of equation-based point and interval forecasts;
* tblCf: output table name of conditionial point and interval forecasts;
* tblCfSim: output table name of simulated conditional forecasts;
proc varmax data=&tblSample.;
model y1 - y3 / p=2 prior noprint;
output alpha=&alpha. lead=&lead. out=&tblF. noprint;
condfore alpha=&alpha. lead=&lead. out=&tblCf. outsim=&tblCfSim.
sdata=&tblScenarios. sid=myscenario
parm=sampling(scenario) nbi=1000 nmc=&nmc.;
run;
%mend;
%macro hcSaveForecasts(dgpi,T,lead,nScenarios,tblDGP,tblF,tblCf,tblAll);
* dgpi: index of DGP;
* T: in-sample sample size;
* lead: future horizons;
* nScenarios: number of scenarios;
* tblDGP: input table name for full-sample data;
* tblF: input table name of equation-based point and interval forecasts;
* tblCf: input table name of conditionial point and interval forecasts;
* tblAll: output table name of point and interval forecasts for all DGPs;
data forecasts;
set &tblDGP.(firstobs=%eval(&T.+1) obs=%eval(&T.+&lead.) keep=Y1 Y2);
* for notation convenience, name equation-based forecasts as S0;
set &tblF.(firstobs=%eval(&T.+1) obs=%eval(&T.+&lead.)
rename=( for1=Y1_S0 for2=Y2_S0
lci1=Y1_LB_S0 uci1=Y1_UB_S0
lci2=Y2_LB_S0 uci2=Y2_UB_S0)
keep=FOR1 LCI1 UCI1 FOR2 LCI2 UCI2);
%do i = 1 %to &nScenarios.;
set &tblCf.(where=(myscenario_S&i.=&i.)
rename=( Y1_MEDIAN=Y1_S&i. Y2_MEDIAN=Y2_S&i.
Y1_LB=Y1_LB_S&i. Y1_UB=Y1_UB_S&i.
Y2_LB=Y2_LB_S&i. Y2_UB=Y2_UB_S&i.
myscenario=myscenario_S&i.)
keep=myscenario Y1_MEDIAN Y2_MEDIAN Y1_LB Y1_UB Y2_LB Y2_UB);
%end;
dgpIndex = &dgpi.;
h = _N_;
drop myscenario_S1 - myscenario_S&nScenarios.;
run;
proc append base=&tblAll. data=forecasts; run;
%mend;
%macro cfEvaluate(tblAll,lead,nSim,nScenarios,qScenario0,scenarioBM,tblEval);
* tblAll: input table name of point and interval forecasts for all DGPs;
* lead: future horizons;
* nSim: number of DGPs;
* nScenarios: number of scenarios;
* qScenario0: whether there is S0 for equation-based forecasts,
1 for yes and 0 for no;
* scenarioBM: the index of the benchmark scenario;
* tblEval: output table name for evaluation results;
proc iml;
use &tblAll.;
read all into d;
close;
lead = &lead.;
nSim = &nSim.;
nScenarios = &nScenarios. + &qScenario0.;
scenarioBM = &scenarioBM.;
nVars = 2;
sMAPE = J(lead,nScenarios*nVars,0);
rsMAPE = J(lead,nScenarios*nVars,0);
sizeCI = J(lead,nScenarios*nVars,0);
rLenCI = J(lead,nScenarios*nVars,0);
do iSim = 1 to nSim;
do h = 1 to lead;
do iScenario = 1 to nScenarios;
do iVar = 1 to nVars;
yF=d[(iSim-1)*lead+h,
nVars+((iScenario-1)*nVars+(iVar-1))*3+1];
yFlb=d[(iSim-1)*lead+h,
nVars+((iScenario-1)*nVars+(iVar-1))*3+2];
yFub=d[(iSim-1)*lead+h,
nVars+((iScenario-1)*nVars+(iVar-1))*3+3];
yFlbBM=d[(iSim-1)*lead+h,
nVars+((scenarioBM-1)*nVars+(iVar-1))*3+2];
yFubBM=d[(iSim-1)*lead+h,
nVars+((scenarioBM-1)*nVars+(iVar-1))*3+3];
y =d[(iSim-1)*lead+h,iVar];
* symmetric Mean Absolute Percentage Error (sMAPE);
if(abs(yF)+abs(y)>0) then do;
sMAPE[h,(iVar-1)*nScenarios+iScenario] =
sMAPE[h,(iVar-1)*nScenarios+iScenario]
+ 2*abs(yF-y)/(abs(yF)+abs(y))/nSim;
end;
* size;
if(y>=yFlb & y<=yFub) then do;
sizeCI[h,(iVar-1)*nScenarios+iScenario] =
sizeCI[h,(iVar-1)*nScenarios+iScenario] + 1/nSim;
end;
* relative length;
if(yFubBM-yFlbBM>0) then do;
rLenCI[h,(iVar-1)*nScenarios+iScenario] =
rLenCI[h,(iVar-1)*nScenarios+iScenario]
+ (yFub-yFlb)/(yFubBM-yFlbBM)/nSim;
end;
end;
end;
end;
end;
do h = 2 to lead;
do iScenario = 1 to nScenarios;
do iVar = 1 to nVars;
sMAPE[h,(iVar-1)*nScenarios+iScenario] =
(sMAPE[h,(iVar-1)*nScenarios+iScenario]
+ sMAPE[h-1,(iVar-1)*nScenarios+iScenario]*(h-1))/h;
end;
end;
end;
do h = 1 to lead;
do iScenario = 1 to nScenarios;
do iVar = 1 to nVars;
* relative symmetric Mean Absolute Percentage Error;
rsMAPE[h,(iVar-1)*nScenarios+iScenario] =
sMAPE[h,(iVar-1)*nScenarios+iScenario]
/ sMAPE[h,(iVar-1)*nScenarios+scenarioBM];
end;
end;
end;
* rearrange results;
n = ncol(sMAPE)/nVars;
evalResult = sMAPE[,1:n];
do iVar = 2 to nVars;
evalResult = evalResult // sMAPE[,(iVar-1)*n+1:iVar*n];
end;
do iVar = 1 to nVars;
evalResult = evalResult // rsMAPE[,(iVar-1)*n+1:iVar*n];
end;
do iVar = 1 to nVars;
evalResult = evalResult // sizeCI[,(iVar-1)*n+1:iVar*n];
end;
do iVar = 1 to nVars;
evalResult = evalResult // rLenCI[,(iVar-1)*n+1:iVar*n];
end;
evalResult = ((1:4)`@J(lead*nVars,1,1))
|| (J(4,1,1)@((1:nVars)`@J(lead,1,1)))
|| (J(4*nVars,1,1)@(1:lead)`)
|| evalResult;
create &tblEval. from evalResult;
append from evalResult;
close;
quit;
%mend;
%macro hcTest(nSim,T,lead,alpha,nmc,nScenarios,qScenario0,scenarioBM,
tblAll,tblEval);
* nSim: number of DGPs;
* T: in-sample sample size;
* lead: future horizons;
* alpha: size of the confidence interval or credible interval;
* nmc: number of Monte Carlo iterations;
* nScenarios: number of scenarios;
* qScenario0: whether there is scenario 0 for equation-based forecasts,
1 for yes and 0 for no;
* scenarioBM: the index of the benchmark scenario;
* tblAll: output table name of point and interval forecasts for all DGPs;
* tblEval: output table name for evaluation results;
%do iSim = 1 %to &nSim.;
%cfSimulateData(&iSim.,&T.,&lead.,t1,t2);
%hcConstructScenarios(&T.,t1,t3);
%hcEstimateAndForecast(t2,t3,&alpha.,&lead.,&nmc.,of,ocf,ocfsim);
%hcSaveForecasts(&iSim.,&T.,&lead.,&nScenarios.,t1,of,ocf,&tblAll.);
%end;
%cfEvaluate(&tblAll.,&lead.,&nSim.,
&nScenarios.,&qScenario0.,&scenarioBM.,
&tblEval.);
%mend;
%let nSim = 1000;
%let T = 200;
%let lead = 4;
%let alpha = 0.05;
%let nmc = 10000;
%let nScenarios = 5;
%let qScenario0 = 1;
%let scenarioBM = 6;
%hcTest(&nSim.,&T.,&lead.,&alpha.,&nmc.,
&nScenarios.,&qScenario0.,&scenarioBM.,
hcForecasts,hcEval);
proc template;
define table hcEvalTemplate;
column col3 col4 col5 col6 col7 col8 col9;
define header hc;
text "Conditional Forecasts, Hard Conditions";
start=col5; end=col8;
end;
define column col3;
header="Horizon";
end;
define column col4;
header="Equation Based"; format=7.5;
end;
define column col5;
header="Scenario 1"; format=12.5;
end;
define column col6;
header="Scenario 2"; format=12.5;
end;
define column col7;
header="Scenario 3"; format=12.5;
end;
define column col8;
header="Scenario 4"; format=12.5;
end;
define column col9;
header="Unconditional"; format=12.5;
end;
end;
run;
%macro cfPrint(template,data);
data _NULL_;
set &data.;
file print ods=(template="&template.");
put _ODS_;
run;
%mend;
%cfPrint(hcEvalTemplate, hcEval(where=(col1=1 and col2=1)));
%cfPrint(hcEvalTemplate, hcEval(where=(col1=1 and col2=2)));
%cfPrint(hcEvalTemplate, hcEval(where=(col1=2 and col2=1)));
%cfPrint(hcEvalTemplate, hcEval(where=(col1=2 and col2=2)));
%cfPrint(hcEvalTemplate, hcEval(where=(col1=3 and col2=1)));
%cfPrint(hcEvalTemplate, hcEval(where=(col1=3 and col2=2)));
%cfPrint(hcEvalTemplate, hcEval(where=(col1=4 and col2=1)));
%cfPrint(hcEvalTemplate, hcEval(where=(col1=4 and col2=2)));
****************************************************************;
* Soft Conditions ;
****************************************************************;
%macro scEstimateAndForecast(tblSample,alpha,lead,nmc,tblUcf,tblUcfSim);
* tblSample: input table name for in-sample data;
* alpha: size of the credible interval;
* lead: future horizons;
* nmc: number of Monte Carlo iterations;
* tblUcf: output table name of unconditional point and interval
forecasts;
* tblUcfSim: output table name of simulated unconditional forecasts;
proc varmax data=&tblSample.;
model y1 - y3 / p=2 prior noprint;
condfore alpha=&alpha. lead=&lead. out=&tblUcf. outsim=&tblUcfSim.
parm=sampling nbi=1000 nmc=&nmc.;
run;
%mend;
%macro scConstructScenarios(T,lead,tblDGP,tblUcf);
* T: in-sample sample size;
* lead: future horizons;
* tblDGP: input table name for full-sample data;
* tblUcf: input table name of unconditional point and interval
forecasts;
* four scenarios are implicitly output:
scenarios i, i=1 to 4: future y3 is known for
lb_j<=y3_j<=ub_j, j=1 to i, where lb_j and ub_j are lower and
upper bounds whose values are saved in the corresponding macro
variables, and y3_j is the j-step-ahead furture value of y3;
data _NULL_;
set &tblDGP.(firstobs=%eval(&T.+1) obs=%eval(&T.+&lead.) keep=Y3);
set &tblUcf.(keep=step Y3_MEDIAN Y3_LB Y3_UB);
array q[&lead.] q1 - q&lead.;
array lb[&lead.] lb1 - lb&lead.;
array ub[&lead.] ub1 - ub&lead.;
retain q lb ub;
if(Y3<=Y3_LB) then do;
q[step] = 1; ub[step] = Y3_LB;
end;
if (Y3>Y3_LB and Y3<=Y3_MEDIAN) then do;
q[step] = 2; lb[step] = Y3_LB; ub[step] = Y3_MEDIAN;
end;
if (Y3>Y3_MEDIAN and Y3<=Y3_UB) then do;
q[step] = 3; lb[step] = Y3_MEDIAN; ub[step] =Y3_UB;
end;
if (Y3>Y3_UB) then do;
q[step] = 4; lb[step] = Y3_UB;
end;
if(_N_=&lead.) then do;
%do i = 1 %to &lead.;
call symputx("lb&i.",lb&i.,'G');
call symputx("ub&i.",ub&i.,'G');
%end;
end;
run;
%mend;
%macro scClassifySimulatedForecasts(lead,tblUcfSim,tblSCSim);
* lead: future horizons;
* tblUcfSim: output table name of simulated unconditional forecasts;
* tblSCSim: output table name of simulated conditional forecasts
under soft conditioins specified in the scenarios;
data &tblSCSim.;
set &tblUcfSim.;
array lb[&lead.] lb1 - lb&lead.;
array ub[&lead.] ub1 - ub&lead.;
array y3f[&lead.] y3_1 - y3_&lead.;
%do j = 1 %to &lead.;
lb[&j.]=&&lb&j.; ub[&j.]=&&ub&j.;
%end;
do myScenario=1 to 4;
outputCond = 1;
do i = 1 to myScenario;
if(outputCond=1) then do;
if(lb[i]=.) then do;
if(y3f[i]<=ub[i]) then outputCond=1;
else outputCond = 0;
end;
else do;
if(ub[i]=.) then do;
if(y3f[i]>lb[i]) then outputCond=1;
else outputCond = 0;
end;
else do;
if(y3f[i]>lb[i] and y3f[i]<=ub[i]) then outputCond=1;
else outputCond = 0;
end;
end;
end;
end;
if(outputCond=1) then output;
end;
run;
proc sort data=&tblSCSim.; by myscenario; run;
%mend;
%macro scGetForecastStats(alpha,lead,tblSCSim,tblSCForecasts);
* alpha: size of the credible interval;
* lead: future horizons;
* tblSCSim: input table name of simulated conditional forecasts
under soft conditioins specified in the scenarios;
* tblSCForecasts: output table name of conditional point and interval
forecasts under soft conditioins specified in the scenarios;
data _NULL_;
lbPctl = &alpha./2*100;
ubPctl = 100-lbPctl;
call symputx("lbPctl",lbPctl,'G');
call symputx("ubPctl",ubPctl,'G');
run;
proc univariate data=&tblSCSim. noprint;
var y1_1 - y1_&lead. y2_1 - y2_&lead.;
output out=oucfx pctlpts=&lbPctl. &ubPctl. 50
pctlpre=%do j=1 %to &lead.; y1_&j. %end;
%do j=1 %to &lead.; y2_&j. %end;
pctlname=_lb _ub _median;
by myscenario;
run;
data &tblSCForecasts.;
set oucfx;
array y1f[&lead.] %do j=1 %to &lead.; Y1_&j._median %end; ;
array y1lb[&lead.] %do j=1 %to &lead.; Y1_&j._lb %end; ;
array y1ub[&lead.] %do j=1 %to &lead.; Y1_&j._ub %end; ;
array y2f[&lead.] %do j=1 %to &lead.; Y2_&j._median %end; ;
array y2lb[&lead.] %do j=1 %to &lead.; Y2_&j._lb %end; ;
array y2ub[&lead.] %do j=1 %to &lead.; Y2_&j._ub %end; ;
do i = 1 to &lead.;
Y1_MEDIAN = y1f[i];
Y1_LB = y1lb[i];
Y1_UB = y1ub[i];
Y2_MEDIAN = y2f[i];
Y2_LB = y2lb[i];
Y2_UB = y2ub[i];
step = i;
output;
end;
keep myscenario step y1_median y1_lb y1_ub y2_median y2_lb y2_ub;
run;
%mend;
%macro scSaveForecasts(dgpi,T,lead,tblDGP,tblUcf,tblSCForecasts,tblAll);
* dgpi: index of DGP;
* T: in-sample sample size;
* lead: future horizons;
* tblDGP: input table name for full-sample data;
* tblUcf: input table name of unconditional point and interval
forecasts;
* tblSCForecasts: input table name of conditional point and interval
forecasts under soft conditioins specified in the scenarios;
* tblAll: output table name of point and interval forecasts for all DGPs;
data forecasts;
set &tblDGP.(firstobs=%eval(&T.+1) obs=%eval(&T.+&lead.) keep=Y1 Y2);
%do i = 1 %to 4;
set &tblSCForecasts.( where=(myscenario_S&i.=&i.)
rename=( Y1_MEDIAN=Y1_S&i. Y2_MEDIAN=Y2_S&i.
Y1_LB=Y1_LB_S&i. Y1_UB=Y1_UB_S&i.
Y2_LB=Y2_LB_S&i. Y2_UB=Y2_UB_S&i.
myscenario=myscenario_S&i.)
keep=myscenario Y1_MEDIAN Y2_MEDIAN Y1_LB Y1_UB Y2_LB Y2_UB);
%end;
set &tblUcf.(
rename=( Y1_MEDIAN=Y1_S5 Y2_MEDIAN=Y2_S5
Y1_LB=Y1_LB_S5 Y1_UB=Y1_UB_S5
Y2_LB=Y2_LB_S5 Y2_UB=Y2_UB_S5 )
keep=Y1_MEDIAN Y2_MEDIAN Y1_LB Y1_UB Y2_LB Y2_UB);
dgpIndex = &dgpi.;
h = _N_;
drop myscenario_S1 - myscenario_S4;
run;
proc append base=&tblAll. data=forecasts; run;
%mend;
%macro scTest(nSim,T,lead,alpha,nmc,nScenarios,qScenario0,scenarioBM,
tblAll,tblEval);
* nSim: number of DGPs;
* T: in-sample sample size;
* lead: future horizons;
* alpha: size of the confidence interval or credible interval;
* nmc: number of Monte Carlo iterations;
* nScenarios: number of scenarios;
* qScenario0: whether there is S0 for equation-based forecasts,
1 for yes and 0 for no;
* scenarioBM: the index of the benchmark scenario;
* tblAll: output table name of point and interval forecasts for all DGPs;
* tblEval: output table name for evaluation results;
%do iSim = 1 %to &nSim.;
%cfSimulateData(&iSim.,&T.,&lead.,t1,t2);
%scEstimateAndForecast(t2,&alpha.,&lead.,&nmc.,oucf,oucfsim);
%scConstructScenarios(&T.,&lead.,t1,oucf);
%scClassifySimulatedForecasts(&lead.,oucfsim,oucfsimx);
%scGetForecastStats(&alpha.,&lead.,oucfsimx,oscf);
%scSaveForecasts(&iSim.,&T.,&lead.,t1,oucf,oscf,&tblAll.);
%end;
%cfEvaluate(&tblAll.,&lead.,&nSim.,
&nScenarios.,&qScenario0.,&scenarioBM.,
&tblEval.);
%mend;
%let nSim = 1000;
%let T = 200;
%let lead = 4;
%let alpha = 0.50;
%let nmc = 100000;
%let nScenarios = 5;
%let qScenario0 = 0;
%let scenarioBM = 5;
%scTest(&nSim., &T.,&lead.,&alpha.,&nmc.,
&nScenarios.,&qScenario0.,&scenarioBM.,
scForecasts,scEval);
proc template;
define table scEvalTemplate;
column col3 col4 col5 col6 col7 col8;
define header sc;
text "Conditional Forecasts, Soft Conditions";
start=col4; end=col7;
end;
define column col3;
header="Horizon";
end;
define column col4;
header="Scenario 1"; format=12.5;
end;
define column col5;
header="Scenario 2"; format=12.5;
end;
define column col6;
header="Scenario 3"; format=12.5;
end;
define column col7;
header="Scenario 4"; format=12.5;
end;
define column col8;
header="Unconditional"; format=12.5;
end;
end;
run;
%cfPrint(scEvalTemplate, scEval(where=(col1=1 and col2=1)));
%cfPrint(scEvalTemplate, scEval(where=(col1=1 and col2=2)));
%cfPrint(scEvalTemplate, scEval(where=(col1=2 and col2=1)));
%cfPrint(scEvalTemplate, scEval(where=(col1=2 and col2=2)));
%cfPrint(scEvalTemplate, scEval(where=(col1=3 and col2=1)));
%cfPrint(scEvalTemplate, scEval(where=(col1=3 and col2=2)));
%cfPrint(scEvalTemplate, scEval(where=(col1=4 and col2=1)));
%cfPrint(scEvalTemplate, scEval(where=(col1=4 and col2=2)));
