The SMC Procedure

This example illustrates the advantages of the fully adapted auxiliary particle filter by using a Gaussian ARCH model.

Consider the following model:

For each distribution, , and .

The fully adapted auxiliary and proposal transition distributions are derived respectively as

First, the decrease in computing time of the APF(ADP) algorithm compared to the APF algorithm is illustrated using the EVALUATE statement. Data sets that have varying numbers of observations are generated. PROC SMC then runs for each data set with varying numbers of particles, and the run time is recorded.

The following statements generate a data set that has a specified number of observations by using the parameter values , , and :

%let lambda_sq = 0.00027;
%let beta0=0.0000495;
%let beta1=0.89275;
%macro generate_data(nobs=);
   proc iml;
      call randseed(123);
      x = j(&nobs,1,.);
      y = j(&nobs,1,.);
      timeID = j(&nobs,1,.);
      x[1]=0;
      call randgen(rny,"Normal");
      y[1]=x[1]+sqrt(&lambda_sq)*rny;
      timeID[1]=1;
      do t=2 to &nobs;
         call randgen(rnx,"Normal");
         x[t]=sqrt(&beta0+&beta1*x[t-1]*x[t-1])*rnx;
         call randgen(rny,"Normal");
         y[t]=x[t]+sqrt(&lambda_sq)*rny;
         timeID[t]=t;
      end;
      create arch var {TimeID x y};
      append ;
      close arch;
   quit;
%mend generate_data;

The following macro repeatedly generates data sets and runs PROC SMC with the fully adapted auxiliary and proposal distributions by using the APF(ADP) and APF algorithms. The average run time and total run time over the 100 replications of the log-likelihood computation are recorded.

%macro runby_npart_nobs(nobs=,N=,alg=,algn=,nb=,);
   %generate_data(nobs=&nobs);
   data mycas.arch;
      set arch;
   run;
   %do r=1 %to 100;
      %let _timer_start = %sysfunc(datetime());
      proc smc data=mycas.arch seed = 123;
         id time = timeID;
         var y;
         statevar x;
         parameters lambda_sq=&lambda_sq beta0=&beta0 beta1=&beta1;
         initstate x ~ degenerate(0) ;
         state x ~ normal(0, sqrt(beta0+beta1*x.l1*x.l1))/proposal=
             normal((beta0+beta1*x.l1*x.l1)*y/(beta0+beta1*x.l1*x.l1+lambda_sq),
             sqrt(lambda_sq*(beta0+beta1*x.l1*x.l1)/
	(beta0+beta1*x.l1*x.l1+lambda_sq)));
         model y ~ normal(x, sqrt(lambda_sq));
         auxmodel y ~ normal(0, sqrt(beta0+beta1*x.l1*x.l1+lambda_sq));
         evaluate nparticle = &N alg=&alg out=mycas.eval;
      run;
      data duration&nb.&r.;
         duration = datetime() - &_timer_start;
      run;
   %end;

   data duration&nb;
      set
      %do r=1 %to 100;
         duration&nb.&r.
      %end;
   ;

   proc univariate data=duration&nb noprint;
      output out= durationsum&algn.&nb. mean=MEAN&algn. sum=SUM&algn.;
   run;

   data durationsum&algn.&nb.;
      set durationsum&algn.&nb;
      N = &N;
      T = &nobs;
   run;
%mend runby_npart_nobs;

The previous macro is then run with , and 1000 and with , and 5000:

%macro runit;
   %let value_nobs = 200 500 1000;
   %let value_npart = 500 1000 2500 5000;
   %local iter iter_nobs iter_npart next_value_nobs next_value_npart;
   %let iter=1; %let iter_nobs=1;  %let iter_npart=1;
   %do %while (%scan(&value_nobs, &iter_nobs) ne );
      %let next_value_nobs = %scan(&value_nobs, &iter_nobs);
      %do %while (%scan(&value_npart, &iter_npart) ne );
         %put  &iter_nobs &iter_npart &iter;
         %let next_value_npart = %scan(&value_npart, &iter_npart);
         %runby_npart_nobs(nobs=&next_value_nobs,N=&next_value_npart,
	alg=%str(APF(ADP)),algn=apf_adp,nb=&iter);
         %runby_npart_nobs(nobs=&next_value_nobs,N=&next_value_npart,
	alg=APF,algn=apf,nb=&iter);
         data durationallsum&iter;
            merge durationsumapf_adp&iter.(in=a) durationsumapf&iter.(in=b);
            by T N;
            if (a=b);
         run;
         %if &iter =1 %then %do;
            %let dataread=durationallsum&iter;
         %end;
         %else %do;
            %let dataread= durationallsum durationallsum&iter;
         %end;
         data durationallsum;
            set &dataread;
         run;
         %let iter_npart = %eval(&iter_npart + 1);
         %let iter = %eval(&iter + 1);
      %end;
      %let iter_nobs = %eval(&iter_nobs + 1);
      %let iter_npart=1;
   %end;
%mend runit;
%runit;

The following statements plot the average log-likelihood computation run time over 100 replications. The results are shown in Figure 12.

proc sgpanel data=durationallsum;
   panelby T / layout=panel columns=1 rows=3 uniscale=column;
   series x=N y=MEANapf/markers markerattrs=(color=blue size=7 symbol=Diamond)
  legendlabel='APF';
   series x=N y=MEANapf_adp/markers markerattrs=(color=red size=7 symbol=asterisk)
  legendlabel='APF(ADP)';
   rowaxis thresholdmin=1 thresholdmax=1 label='Time (seconds)';
   colaxis label='Number of Particles'
        valueshint values=(500 1000 2500 5000);
run;

The following statements plot the total log-likelihood computation run time over 100 replications by using a single computing node. The results are shown in Figure 13.

proc sgpanel data=durationallsum;
   panelby T / layout=panel columns=1 rows=3 uniscale=column;
   series x=N y=SUMapf/markers markerattrs=(color=blue size=7 symbol=Diamond)
  legendlabel='APF';
   series x=N y=SUMapf_adp/markers markerattrs=(color=red size=7 symbol=asterisk)
  legendlabel='APF(ADP)';
   rowaxis thresholdmin=1 thresholdmax=1 label='Time (seconds)';
   colaxis label='Number of Particles'
        valueshint values=(500 1000 2500 5000);
run;

Figure 12 and Figure 13 show that the APF(ADP) algorithm is faster than the APF algorithm with fully adapted auxiliary and proposal distributions, especially when the number of particles and/or the number of time points is large. A similar conclusion can be made about the FILTER, FORECAST, LEARN, and SMOOTH statements.

In addition to being faster, the APF(ADP) algorithm results in a smaller variance for the estimated log likelihood and a smaller optimal number of particles compared to APF and BF. The next code aims to illustrate this.

The following statements simulate 500 observations on the basis of the defined ARCH model and load the data set arch into a data table in your CAS session that is associated with the mycas CAS engine libref:

%let lambda_sq = 0.00027;
%let beta0=0.0000495;
%let beta1=0.89275;
proc iml;
   call randseed(123);
   x = j(500,1,.);
   y = j(500,1,.);
   timeID = j(500,1,.);
   x[1]=1;
   call randgen(rny,"Normal");
   y[1]=x[1]+sqrt(&lambda_sq)*rny;
   timeID[1]=1;
   do t=2 to 500;
     call randgen(rnx,"Normal");
     x[t]=sqrt(&beta0+&beta1*x[t-1]*x[t-1])*rnx;
     call randgen(rny,"Normal");
     y[t]=x[t]+sqrt(&lambda_sq)*rny;
     timeID[t]=t;
   end;
   create arch var {TimeID x y};
   append ;
   close arch;
quit;

data mycas.arch;
   set arch;
run;

Following the guideline for choosing the optimal number of particles provided by Pitt et al. (2012), the next group of statements estimate the log likelihood 100 times for a given number of particles and filtering algorithm, and then they compute the standard deviation of the log likelihood over the 100 replications:

%macro repeat_run(npart=,alg=,algn=,nb=);
   %do i=1 %to 100;
      proc smc data=mycas.arch seed = &i;
         id time = timeID;
         var y;
         statevar x;
         parameters lambda_sq=&lambda_sq beta0=&beta0 beta1=&beta1;
         initstate x ~ degenerate(1) ;
         state x ~ normal(0, sqrt(beta0+beta1*x.l1*x.l1))
	%if %length(&alg) gt 3 %then /proposal=
	normal((beta0+beta1*x.l1*x.l1)*y/(beta0+beta1*x.l1*x.l1+lambda_sq),
	sqrt(lambda_sq*(beta0+beta1*x.l1*x.l1)/(beta0+beta1*x.l1*x.l1+lambda_sq)));;
         model y ~ normal(x, sqrt(lambda_sq));
         auxmodel y ~ normal(0, sqrt(beta0+beta1*x.l1*x.l1+lambda_sq));
         evaluate nparticle = &npart alg=&alg out=mycas.eval&nb.&i.;
      run;
      data eval&nb.&i.(keep = LogLikelihood);
         set mycas.eval&nb.&i. end=eof;
         if eof=1;
      run;
   %end;
   data evalall&nb;
      set
      %do i=1 %to 100;
         eval&nb.&i.
      %end;
   ;

   proc univariate data= evalall&nb vardef=n noprint;
      output out= evalallsum&algn.&nb.  mean=MEAN std=STDDEV var=VARIANCE n=R;
   run;

   data evalallsum&algn.&nb.;
      set evalallsum&algn.&nb(drop=R);
      N = ∂̸
      algorithm = "&alg";
   run;
%mend repeat_run;

The following statements run the macro with , and 500 particles for the APF(ADP), APF, and BF algorithms:

%macro runit2;
   %let value = 15 20 25 50 75 90 100 200 300 400 500;
   %local iter next_value;
   %let iter=1;
   %do %while (%scan(&value, &iter) ne );
      %let next_value = %scan(&value, &iter);
      %repeat_run(npart=&next_value,alg=%str(APF(ADP)),algn=1,nb=&iter);
      %repeat_run(npart=&next_value,alg=APF,algn=2,nb=&iter);
      %repeat_run(npart=&next_value,alg=BF,algn=3,nb=&iter);
      data evalallsum&iter;
         set evalallsum1&iter evalallsum2&iter evalallsum3&iter;
      run;
      %if &iter =1 %then %do;
         %let dataread=evalallsum&iter;
      %end;
      %else %do;
         %let dataread= evalallsum evalallsum&iter;
      %end;

      data evalallsum;
         set &dataread;
      run;
      %let iter = %eval(&iter + 1);
   %end;
%mend runit2;
%runit2;

The log-likelihood standard deviation across the different number of particles and filtering algorithms is then plotted using the following code. The results are shown in Figure 14. Following Pitt et al. (2012), a reference line corresponding to a standard deviation of 0.92 is added to help determine the optimal number of particles.

data evalallsum;
   set evalallsum;
   ref = 0.92;
run;

proc sort data=evalallsum;
   by algorithm N;
run;
proc sgpanel data=evalallsum;
   title ;
   panelby algorithm / layout=panel columns=3 rows=1 uniscale=all;
   series x=N y=STDDEV/markers markerattrs=(color=blue size=7 symbol=plus)
     legendlabel='StdDev ';
   series x=N y=ref / legendlabel='Cutoff';
   rowaxis thresholdmin=1 thresholdmax=1 label='Standard Deviation' max=25;
   colaxis label='Number of Particles'
        valueshint
        values=(0 15 20 25 50 75 90 100 200 300 400 500);
   keylegend/position=bottom;
run;

Figure 14 shows that the variability in the log likelihood is significantly smaller for APF(ADP) than for APF and BF. Moreover, the optimal numbers of particles that are determined as the number of particles corresponding to a standard deviation of 0.92 (intersection of the red reference line and standard deviation line) are 20, 450, and 450 for APF(ADP), APF, and BF, respectively. Because ADF(ADP) requires fewer particles in order to reach the optimal standard deviation, this further improves the speed of computation.

A final advantage of APF(ADP) over APF and BF is its lower Markov chain autocorrelation. To illustrate this, particle Markov chain Monte Carlo (MCMC) is run using the arch data set and the APF(ADP), APF, and BF algorithms.

The parameters and are transformed so that the new parameters and are unconstrained: and . The following prior distributions are specified:

The following macro runs particle MCMC and saves the acceptance rate and diagnostics summary tables for a given filtering algorithm and number of particles:

%macro run_learn(npart=,alg=,algn=);
   ods output DiagSummaries= DiagSummary_&algn.
         AccRates=AccRates_&algn.;
   proc smc data=mycas.arch seed =123;
      id time = timeID;
      var y;
      statevar x;
      parameters lambda_sq=0.005 theta0=-6 theta1=-0.3;
      initstate x ~ degenerate(1) ;
      state x ~ normal(0, sqrt(exp(theta0)+exp(theta1)*x.l1*x.l1))
%if %length(&alg) gt 3 %then /proposal=
         normal((exp(theta0)+exp(theta1)*x.l1*x.l1)*y/(exp(theta0)
	+exp(theta1)*x.l1*x.l1+lambda_sq),
         sqrt(lambda_sq*(exp(theta0)+exp(theta1)*x.l1*x.l1)/(exp(theta0)+
	exp(theta1)*x.l1*x.l1+lambda_sq)));;
      model y ~ normal(x, sqrt(lambda_sq));
      auxmodel y ~ normal(0, sqrt(exp(theta0)+exp(theta1)*x.l1*x.l1+lambda_sq));
      prior lambda_sq ~ igamma(0.01,0.01);
      prior theta0 ~ normal(-10,4);
      prior theta1 ~ normal(-0.13,2.8);
      learn nparticle =  &npart algorithm=&alg
         method = PMCMC(nsample = 50000 nbi=10000 samsum samset seed = 1
         thin=1 outpost = mycas.outle_&algn.
         statistics=(corr cov interval(alpha = (0.01))
         summary(perc = (0.5 2.5 25 50 75 97.5 99)))
         diagnostics(settings)=all sampler = rwm(ntu=700 ntunestage=48)) ;
   run;

   data AccRates_&algn.;
      set AccRates_&algn.;
      Algorithm = "&alg";
      N = ∂̸
   run;
%mend run_learn;

The macro is then run for the APF(ADP), APF, and BF algorithms by using the optimal number of particles determined previously. A computing time variable is defined as .

%run_learn(npart=20,alg=%str(APF(ADP)),algn=apf_adp);
%run_learn(npart=450,alg=APF,algn=apf);
%run_learn(npart=450,alg=BF,algn=bf);

data AccRates;
   set AccRates_apf_adp AccRates_bf AccRates_apf;
run;

data DiagSummary_bf;
   set DiagSummary_bf;
   ComputingTime = 450*CorrTime/1000;
run;
proc transpose data=DiagSummary_bf out=DiagSummary_bf_w;
   id Parameter;
run;
data DiagSummary_bf_w; set DiagSummary_bf_w;
   order=2; Algorithm='BF'; N = 450;run;

data DiagSummary_apf;
   set DiagSummary_apf;
   ComputingTime = 450*CorrTime/1000;
run;
proc transpose data=DiagSummary_apf out=DiagSummary_apf_w;
   id Parameter;
run;
data DiagSummary_apf_w; set DiagSummary_apf_w;
   order=1; Algorithm='APF'; N = 450;run;

data DiagSummary_fapf;
   set DiagSummary_fapf;
   ComputingTime = 20*CorrTime/1000;
run;
proc transpose data=DiagSummary_fapf out=DiagSummary_fapf_w;
   id Parameter;
run;
data DiagSummary_fapf_w; set DiagSummary_fapf_w;
   order=3; Algorithm='APF(ADP)'; N = 20;run;

data DiagSummary;
   set DiagSummary_fapf_w DiagSummary_bf_w DiagSummary_apf_w;
run;

The following statements print the acceptance rates for each algorithm, as shown in Output 23.5.1. The acceptance rates of the three algorithms are around the 0.234 optimal value for the random walk Metropolis sampler, despite APF(ADP) using far fewer particles.

proc sort data=AccRates; by descending N;run;
proc print data=AccRates noobs;
run;

The following statements print the autocorrelation time and the indicator variable that demonstrates whether autocorrelation remains in the chain, as shown in Output 23.5.2 and Output 23.5.3, respectively:

data CorrTime;
   set DiagSummary;
   where _NAME_ = 'CorrTime';
run;
proc sort data=CorrTime; by order;run;
proc print data=CorrTime noobs;
   var Algorithm N lambda_sq theta0 theta1;
run;
data AutoCorrRemains;
   set DiagSummary;
   where _NAME_ = 'AutoCorrRemains';
run;
proc sort data=AutoCorrRemains; by order;run;
proc print data=AutoCorrRemains noobs;
   var Algorithm N lambda_sq theta0 theta1;
run;

Output 23.5.2 and Output 23.5.3 show that significant autocorrelation still remains in the Markov chain for and when the APF and BF algorithms are used, indicating that you need to run the Markov chain longer or increase the tuning time of the proposal distribution.

Finally, the following statements print the computing time, which is shown in Output 23.5.4:

data ComputingTime;
   set DiagSummary;
   where _NAME_ = 'ComputingTime';
run;
proc sort data=ComputingTime; by order;run;
proc print data=ComputingTime noobs;
   var Algorithm N lambda_sq theta0 theta1;
run;

Output 23.5.4 shows that APF(ADP) has a lower computing time, defined as , than APF and BF for all the parameters in the model. Clearly, the APF(ADP) algorithm with 20 particles is more efficient than the APF or BF algorithm with 450 particles.