DYNAMICLINEAR Procedure

This example shows how the DYNAMICLINEAR procedure enables you to track volatility by using the simultaneous graphical dynamic linear model (SGDLM) and how the parental set of the SGDLM can help you predict volatility more accurately.

Using the same notation as in the section Simultaneous Graphical Dynamic Linear Models, consider generating m-dimensional time series that follow the SGDLM by using the following steps:

This example assumes that each variable has five elements in its parental set. The following statements simulate 300 observations () for 11 time series () with the volatility factors and :

*macro variables;
%let seed = 12345;
%let m = 11;
%let T = 300;
%let nparent = 5;
%let sigma_phi = 0.0001;
%let sigma_gamma = 0.01;

proc iml;
   call randseed(&seed.);
   *define the parental sets;
   nparent_j =j(&m,1,0);
   do h=1 to &m.;
      nparent_j[h] = &nparent.;
   end;

   parents=j(&m.,&m.,0);
   all_ind = 1:&m.;
   do h=1 to &m.;
      p_gam=nparent_j[h];
      if p_gam >0 then do;
         all_ind_j = all_ind[loc(all_ind ^= h)];
         parent_ju = sample(all_ind_j, p_gam, "NoReplace");
         parent_ju_c = parent_ju`;
         call sort(parent_ju_c);
         parent_j = parent_ju_c`;
         *print p_gam parent_j;
         parents[h,parent_j]=1;
      end;
   end;

   parentsc= {};
   all_ind = 1:&m.;
   do h=1 to &m.;
      parents_h = loc(parents[h,]=1);
      npar_h = ncol(parents_h);
      if npar_h >0 then do;
         do k=1 to npar_h;
            pjk = h||parents_h[k];
            parentsc = parentsc//pjk;
         end;
      end;
   end;

   *generate data;
   y=j(&T,&m.+1);
   do t=1 to &T;
      gamma_t = j(&m, &m, 0);
      all_ind = 1:&m.;
      do h=1 to &m.;
         p_gam=nparent_j[h];
         if p_gam >0 then do;
            parent_j = all_ind[loc(parents[h,] = 1)];
            gamj = j(p_gam,1);
            call randgen(gamj,"normal");
            gamma_t[h,parent_j`]=&sigma_gamma.*gamj`;
         end;
      end;

      *generate mut;
      mu_ts = j(&m, 1);
      call randgen(mu_ts, "Normal");
      mu_t = &sigma_phi.*mu_ts;

      *generate Lambdat;
      lambda_t = j(&m, 1);
      call randgen(lambda_t, "Gamma",1);
      lambda_t = inv(diag(lambda_t));

      *compute A_t, Sigma_t;
      A_t = inv(I(&m)-gamma_t);
      mean_t = A_t*mu_t;
      Sigma_t = A_t*inv(lambda_t)*A_t`;
      /* the Cholesky root of the full covariance matrix */
L_t = root(Sigma_t);
      z_t = j(&m,1);
      call randgen(z_t, "Normal"); /* z is MVN(0,I(&m)) */
      y_t = mean_t + L_t`*z_t;
      y[t,1:&m.]=y_t`;
      y[t,&m.+1]=t;
   end;

   yr = y[loc(y[,&m.+1] <= &T-1),];

   yvarNames = "y1":"y&m.";
   varNames =yvarNames||"t";
   create data_correlated from y [colname=varNames];
   append from y;
   close data_correlated;

   create data_correlatedr from yr [colname=varNames];
   append from yr;
   close data_correlatedr;

   create parents from parents [colname=yvarNames];
   append from parents;
   close parents;

   varNamesc = {'ParentID', 'ChildID'};
   create adj_table from parentsc [colname=varNamesc];
   append from parentsc;
   close adj_table;
run;

proc sort data=data_correlatedr;
   by t;
run;

The following statements load the generated data set data_correlatedr and the defined adjacency matrix data set adj_table into data tables in your session that is associated with the mylib libref. The DATA step assumes that your libref is named mylib, but you can substitute any appropriately defined libref.

data mylib.test_data;
   set data_correlatedr;
run;

data mylib.adj_table;
   set adj_table;
run;

Initialization settings are very important in SGDLMs. PROC DYNAMICLINEAR uses the identity matrix as the default initialization covariance. However, the identity matrix does not work well in this example. Following Gruber and West (2016), the covariance matrix is chosen as the initialization covariance matrix for all variables in the SGDLM. The following DATA step inputs the covariance matrices into the data set initial_cov_with_adj and then loads it into a data table in your session that is associated with the mylib libref:

data initial_cov_with_adj;
   input VarID ParmID ParmValue;
datalines;
1 1 0.0001
1 2 0.01
1 3 0.01

   ... more lines ...   

11 3 0.01
11 4 0.01
11 5 0.01
11 6 0.01
;

data mylib.initial_cov_with_adj;
   set initial_cov_with_adj;
run;

The following statements use PROC DYNAMICLINEAR to define the preceding SGDLM and then generate predictions for the simulated data set:

proc dynamicLinear data=mylib.test_data nSimulations=10000 seed=12345;
   id t;
   model y1-y11;
   parentalSet inAdjacency=mylib.adj_table;
   filter out(percentiles=(0.5 2.5 50 97.5 99.5))=mylib.out_filter_with_adj
      outFilterInfo=(klDivergence=mylib.kldiv)
      back=0;
   forecast out(percentiles=(0.5 2.5 50 97.5 99.5))=mylib.out_forecast_with_adj;
   initialDistribution inCovariance(diagonal)=mylib.initial_cov_with_adj
      shape=(default=5) rate=(default=1.5);
   discountFactor beta=(default=0.98) deltaPhi=(default=0.98)
      deltaGamma=(default=0.99);
run;

The summary tables of the model are shown in Output 17.1.1.

SGDLMs use the Kullback-Leibler (KL) divergence to measure the difference between the estimated and true posterior distributions. The KL divergence is calculated during the importance sampling of posterior updating. PROC DYNAMICLINEAR provides the KLDIVERGENCE= suboption of the OUTFILTERINFO= option in the FILTER statement to output the values of the KL divergence and the effective sample size at each step. The following statements create a plot of the two values at each step. The plot is shown in Output 17.1.2.

data kldiv;
   set mylib.kldiv;
run;

proc sort data=kldiv;
   by t;
run;

proc sgplot data=kldiv;
   series x=t y=KLDivergence / legendlabel='KL Divergence Values'
      lineattrs=(color=blue pattern=solid);
   series x=t y=EffectiveSampleSize / legendlabel='Effective Sample Size Values'
      lineattrs=(color=red pattern=shortdash) Y2Axis;
   label t="Step" KLDivergence="KL Divergence"
      EffectiveSampleSize="Effective Sample Size";
run;

Output 17.1.2 shows that PROC DYNAMICLINEAR creates a good estimation of the posterior distribution in this example: first, the effective sample size converges to the number of simulations (10,000); then the KL divergence converges to 0. For more information about the quality of the posterior estimation in SGDLMs, see Gruber and West (2016).

The following statements run the same SGDLM without the parental sets to compare the results to the previous ones. The purpose of this comparison is to demonstrate the effect of the parental sets. The statements first define the initial covariance matrices and then run PROC DYNAMICLINEAR. All other settings are the same as those in the preceding statements.

data mylib.initial_cov_wo_adj;
   input VarID ParmID ParmValue;
datalines;
1 1 0.0001
2 1 0.0001
3 1 0.0001
4 1 0.0001
5 1 0.0001
6 1 0.0001
7 1 0.0001
8 1 0.0001
9 1 0.0001
10 1 0.0001
11 1 0.0001
;

proc dynamicLinear data=mylib.test_data;
   id t;
   model y1-y11;
   filter out(percentiles=(0.5 2.5 50 97.5 99.5))=mylib.out_filter_wo_adj
   back=0;
   forecast out(percentiles=(0.5 2.5 50 97.5 99.5))=mylib.out_forecast_wo_adj;
   initialDistribution inCovariance(diagonal)=mylib.initial_cov_wo_adj
      shape=(default=5) rate=(default=1.5);
   discountFactor beta=(default=0.98) deltaPhi=(default=0.98)
      deltaGamma=(default=0.99);
run;

The following statements create a plot of the 60-step moving average values, which are the mean predictions with and without the parental sets of the variable y1 at each step, starting from . The moving average values of the variable y1 are calculated by the EXPAND procedure and recorded in the data table test_data_ma. Then the data are merged and recorded in the data table y1_mean_comb. The plot is shown in Output 17.1.3.

proc expand data=data_correlatedr out=test_data_ma method=none;
   id t;
   convert y1 = y1_ma / transout=(movave 60);
run;

data y1_mean_comb;
   format t y1_ma y1_adj y1_wo_adj Best12.;
   set test_data_ma(keep=t y1_ma where=(t>=60));
   set mylib.out_filter_with_adj(keep=t Mean_y1 where=(t>=59)
      rename=(Mean_y1=y1_adj));
   set mylib.out_filter_wo_adj(keep=t Mean_y1 where=(t>=59)
      rename=(Mean_y1=y1_wo_adj));
   t = t+1;
run;

proc sort data=y1_mean_comb;
   by t;
run;

proc sgplot data=y1_mean_comb;
   series x=t y=y1_ma / legendlabel='Moving Average Values of y1'
      markers markerattrs=(color=red size=5 symbol=circlefilled)
      lineattrs=(color=red pattern=shortdash);
   series x=t y=y1_adj / legendlabel='Mean_y1 with the Parental Sets'
      markers markerattrs=(color=blue size=5 symbol=triangle)
      lineattrs=(color=blue pattern=solid);
   series x=t y=y1_wo_adj / legendlabel='Mean_y1 without the Parental Sets'
      markers markerattrs=(color=cyan size=5 symbol=circle)
      lineattrs=(color=cyan pattern=solid);
   label t="Step" y1_ma="Values of the Means";
run;

Output 17.1.3 shows that in this example, the SGDLM with the parental sets and the SGDLM without the parental sets provide similar mean predictions.

The following statements create a plot of the 60-step moving average variances, which are the variance predictions with and without the parental sets of the variable y1 at each step, starting from . The moving average variances of the variable y1 are calculated by the EXPAND procedure and recorded in the data table test_data_ma_var. Then the data are merged and recorded in the data table y1_var_comb. The plot is shown in Output 17.1.4.

/*calculate 60-step moving average variances of the variable y1*/
proc expand data=data_correlatedr out=test_data_ma_var method=none;
   id t;
   convert y1 = y1_ma_var / transout=(MOVVAR 60);
run;

data y1_var_comb;
   format t y1_ma_var y1_adj y1_wo_adj Best12.;
set test_data_ma_var(keep=t y1_ma_var where=(t>=60));
set mylib.out_filter_with_adj(keep=t Variance_y1 where=(t>=59)
      rename=(Variance_y1=y1_adj));
set mylib.out_filter_wo_adj(keep=t Variance_y1 where=(t>=59)
      rename=(Variance_y1=y1_wo_adj));
   t = t+1;
run;

proc sort data=y1_var_comb;
   by t;
run;

proc sgplot data=y1_var_comb;
   series x=t y=y1_ma_var / legendlabel='Moving Average Variances of y1'
      markers markerattrs=(color=red size=5 symbol=circlefilled)
      lineattrs=(color=red pattern=shortdash);
   series x=t y=y1_adj / legendlabel='Variance_y1 with the Parental Sets'
      markers markerattrs=(color=blue size=5 symbol=triangle)
      lineattrs=(color=blue pattern=solid);
   series x=t y=y1_wo_adj / legendlabel='Variance_y1 without the Parental Sets'
      markers markerattrs=(color=cyan size=5 symbol=circle)
      lineattrs=(color=cyan pattern=solid);
   label t="Step" y1_ma_var="Values of the Variances";
run;

Output 17.1.4 shows that the SGDLM with the parental sets provides better volatility tracking than the SGDLM without the parental sets. The results are consistent with the conclusions of Gruber and West (2016).