The SMC Procedure

This example illustrates how to define an autoregressive stochastic volatility model in the SMC procedure and then estimate the parameters by using the particle Markov chain Monte Carlo (PMCMC) method.

The stochastic volatility model assumes that the volatility is a latent variable with unexpected innovation terms (Shephard 2005). The discrete-time autoregressive stochastic volatility (ARSV) model proposed by Taylor (1982) is one of the simplest stochastic volatility models. The ARSV model, in which the log variance follows a linear Gaussian process, is given as

where denotes the log daily return series and denotes the daily volatility series. Suppose that the true parameter values are , , , and . In order to bring these initial settings in line with the real data, here the denotes the 100-times-magnified log daily return series; that is, , where is the asset price at day t ().

The following statements simulate 280 observations by using the given stochastic volatility model:

/*--- parameter values ---*/
%let phi = 0.9;
%let gamma = 0.4;
%let beta = 1;

/*--- Simulate a sample for the stochastic volatility model ---*/
data asv_all;
   xlag = 0.45;
   do t = 1 to 280;
      eps = normal(1);
      x = &phi*xlag + &gamma*eps;
      eta = normal(1);
      y = &beta*exp(x/2)*eta;
      output;
      xlag = x;
   end;
run;

The following DATA steps load the data set asv_all into a data table in your CAS session that is associated with the mycas CAS engine libref. The DATA step assumes that your CAS engine libref is named mycas, but you can substitute any appropriately defined CAS engine libref. The data set asv_all is split such that the first 275 observations are used as the training set and the last 5 observations are used as the test set:

data asv_train;
   set asv_all;
   if (t le 275) then output;
run;

data asv_test;
   set asv_all;
   if (t > 275) then output;
run;

data mycas.asv_train;
   set asv_train;
run;

The model specification defines the transition and measurement distributions. For the initialization distribution of the state variable, the normal distribution whose mean is 0 and variance is 4 is used. Moreover, the prior distributions for the parameters are defined as

The following statements use PROC SMC to define the preceding ARSV model and then use the PMCMC method to estimate the three parameters for the simulated data set:

proc smc data=mycas.asv_train seed=1234;
   id time = t;
   var y;
   statevar x;
   prior phi~uniform(0,1);
   prior gamma_sq beta_sq~igamma(0.01,0.01);
   parm phi=0.6 gamma_sq=0.9 beta_sq=1.5;
   initstate x~normal(0,2);
   state x~normal(x.l1*phi, sqrt(gamma_sq));
   model y~normal(0, sqrt(beta_sq)*exp(x/2));
   auxmodel y~normal(0, sqrt(beta_sq)*exp(phi*x.l1/2));
   learn nparticle=256 algorithm=apf method=PMCMC(nsample=25000
      nchain=8 nbi=5000 thin=4 seed=1
      outpost=mycas.out samset samsum
      sampler=rwm(ntu=250)
      statistics=all diagnostics=all);
   score outmodel=mycas.asvmodel;
run;

For the three parameters, is represented by phi, is represented by gamma_sq, and is represented by beta_sq. A state variable, x, an a dependent (observation) variable, y, are also specified according to the model specification.

The model information and parameter learning results are stored in the mycas.asvmodel. As for the PMCMC, the likelihood is estimated by the particle filter and the prior density is specified in the PRIOR statement.

The summary tables of the model information are shown in Figure 1.

The summary tables of the posterior samples are shown in Figure 2.

The correlation and covariance matrices of the parameter samples that are generated by the PMCMC method are shown in Figure 3.

The autocorrelation statistics, Heidelberger-Welch diagnostics, effective sample size (ESS) diagnostics, and Monte Carlo Standard Error (MCSE) diagnostics for each parameter are shown in Figure 4.

The following statements draw the histogram of the posterior samples of phi, which is shown in Figure 5:

proc sgplot data=mycas.out;
   histogram phi;
   refline 0.9 /axis=x lineattrs=(thickness=5 color=red) label="True Value";
run;

The following statements draw the histogram of the posterior samples of gamma_sq, which is shown in Figure 6:

proc sgplot data=mycas.out;
   histogram gamma_sq;
   refline 0.16 /axis=x lineattrs=(thickness=5 color=red) label="True Value";
run;

The following statements draw the histogram of the posterior samples of beta_sq, which is shown in Figure 7:

proc sgplot data=mycas.out;
   histogram beta_sq;
   refline 1 /axis=x lineattrs=(thickness=5 color=red) label="True Value";
run;

You can apply the model information and parameter estimates from the data table mycas.asvmodel to the original data set by specifying the INMODEL= option in the SCORE statement. The posterior mean acquired by the PMCMC method are used as the parameters in the smoothing and forecasting process. The following statements generate the smoothing and forecasting estimates:

proc smc data=mycas.asv_train seed=1 nthreads=32;
   smooth nparticle=10240 parm=median alg=sir method=fixedlag(25)
      out(pct=2.5,97.5)=mycas.asvsm;
   forecast nparticle=10240 parm=median alg=sir lead=5
      out(pct=2.5,97.5)=mycas.asvfc;
   score inmodel=mycas.asvmodel;
run;

The following statements sort the summary table of the smoothing estimates of the state variable, x, in the CAS session and create a relevant data table, smoothRes, in the Work library to record the smoothing estimates and actual values:

data asvsm; set mycas.asvsm(where=(VariableID=1 and TimeID>=1)); run;
proc sort data=asvsm; by timeId; run;

data smoothRes;
   set asvsm (keep=timeId Q2 P_1 P_2);
   set asv_train (keep=x t);
run;

The following statements draw the plot of the smoothing median estimates and the actual values of x. The output is shown in Figure 8, where the shadow band indicates the range between the 2.5th and 97.5th percentiles of the smoothing estimates.

proc sgplot data=smoothRes;
   band x=TimeID upper=P_2 lower=P_1 / transparency=.3 legendlabel='95% CI';
   series x=TimeId y=Q2 /legendlabel='Smoothing Median'
      markers markerattrs=(color=blue size=5 symbol=triangle)
      lineattrs=(color=blue pattern=solid);
   series x=TimeId y=x /legendlabel='Actual Value'
      markers markerattrs=(color=red size=5 symbol=circlefilled)
      lineattrs=(color=red pattern=shortdash);
   yaxis label='x, Smoothing';
   xaxis label='Time';
run;

The following statements sort the summary table of the forecasting estimates of x in the CAS session and create a relevant data table, forecastRes_x, in the Work library to record the forecasting estimates and actual values in the forecasting stage:

data asvfc_x; set mycas.asvfc(where=(VariableID=1)); run;
proc sort data=asvfc_x; by _t; run;

data forecastRes_x;
   set asvfc_x (keep=Q2 P_1 P_2 _t);
   set asv_test (keep=x t);
run;

The following statements draw the plot of the forecasting median estimates and the actual values of x. The output is shown in Figure 9, where the shadow band indicates the range between the 2.5th and 97.5th percentiles of the forecasting estimates.

proc sgplot data=forecastRes_x;
   band x=t upper=P_2 lower=P_1 /transparency=.3 legendlabel='95% CI';
   series x=t y=Q2 /legendlabel='Forecasting Median'
      markers markerattrs=(color=blue size=5 symbol=triangle)
      lineattrs=(color=blue pattern=solid);
   scatter x=t y=x /legendlabel='Actual Value'
      markerattrs=(color=red size=5 symbol=circlefilled);
   xaxis integer label='Time';
   yaxis label='x, Forecasting';
run;

The following statements sort the summary table of the forecasting estimates of the dependent variable, y, in the CAS session and create a relevant data table, forecastRes_y, in the Work library to record the forecasting estimates and actual values in the forecasting stage:

data asvfc_y; set mycas.asvfc(where=(VariableID=2)); run;
proc sort data=asvfc_y; by _t; run;

data forecastRes_y;
   set asvfc_y (keep=Q2 P_1 P_2 _t);
   set asv_test (keep=y t);
run;

The following statements draw the plot of the forecasting median estimates and the actual values of y. The output is shown in Figure 10, where the shadow band indicates the range between the 2.5th and 97.5th percentiles of the forecasting estimates.

proc sgplot data=forecastRes_y;
   band x=t upper=P_2 lower=P_1 /transparency=.3 legendlabel='95% CI';
   series x=t y=Q2 /legendlabel='Forecasting Median'
      markers markerattrs=(color=blue size=5 symbol=triangle)
      lineattrs=(color=blue pattern=solid);
   scatter x=t y=y /legendlabel='Actual Value'
      markerattrs=(color=red size=5 symbol=circlefilled);
   xaxis integer label='Time';
   yaxis label='y, Forecasting';
run;