The SMC Procedure

This example illustrates how to define a stochastic susceptible-infected-recovered (SSIR) model with a constant transmission rate in the SMC procedure and then estimate the parameters using the PMCMC method.

Consider the following SSIR model within the framework of the state space model:

Here , a Poisson distributed instance, denotes the number of new infections at time ; denotes the number of susceptible individuals at time t; is the number of new observed infections at time t; is the cumulative number of recovered infections until time t; is the number of documented active infectious individuals at time t; and is the number of undocumented active infectious individuals at time t. Moreover, N denotes the population, r denotes the basic reproduction number, p is the diagnosis rate, D denotes the average time to removal after infection, denotes the recovery rate, and is the transmission rate. In this SSIR model, there are three parameters: r, D, and p; and there are five state variables: , S, , , and .

The prior specification for the parameters are given as follows:

The initial values and bounds for the parameters are given by

where the bounds of r and p are identical to the supports of their prior distributions, whereas the bounds of D are different from the support of its prior distribution, which is . To achieve this, the bounds can be specified in the PARAMETER statement as a supplement to the supports of the prior distributions.

The following statements simulate 89 observations based on the given stochastic volatility model. The first 75 observations are used as a training data set, and the remaining 14 observations are used as a test data set.

*macro variables;
%let N=20000000; *population;
%let Y_1 = 100;
%let IU_1 = 500;
%let RC_1 = 0;
%let D = 8.7;
%let r = 3.1;
%let p = 0.25;
data sir_all;
   call streaminit(1);
   t = 1;
   ID = &Y_1;
   IU = &IU_1;
   RC = &RC_1;
   S = &N - ID - IU;
   xi = RAND('POISSON', &r/(&N*&D)*S*(ID+IU));
   Y = &Y_1;
   sum = S + IU + ID + RC;
   output;
   do t = 2 to 89;
      RC = RC + 1/&D*(IU + ID);
      IU = (1-1/&D)*IU + xi - Y;
      ID = (1-1/&D)*ID + Y;
      S = S - xi;
      xi = RAND('POISSON', &r/(&N*&D)*S*(ID+IU));
      Y = RAND('BINOMIAL', &p, int(IU));
      sum = S + IU + ID + RC; *sum is always equal to N;
      output;
   end;
run;

The following DATA steps load the first 75 observations of the data set sir_all into a data set sir_train and then load that data set 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 last 14 observations of the data set sir_all are loaded into a data set sir_test.

data sir_train;
   set sir_all;
   if (t le 75) then output;
run;

data sir_test;
   set sir_all;
   if (t > 75) then output;
run;

data mycas.sir_train;
   set sir_train;
run;

The following statements use PROC SMC to define the preceding SSIR model and then use the PMCMC method to estimate the three parameters for the simulated data set. The posterior medians will be used as the parameters in the smoothing and forecasting process.

proc smc data=mycas.sir_train seed = 1234 nthreads = 32;
   id time = t;
   var y;
   statevar S ID IU xi RC;
   smooth nparticle=20480 parm=median alg=sir method=fixedlag(16)
      out(pct=2.5,97.5)=mycas.sm;
   forecast nparticle=20480 parm=median alg=sir lead=14
      out(pct=2.5,97.5)=mycas.fc;
   parm r=2.9, D=9.3, p=0.3 /lower=(2.6 6.5 0.1) upper=(3.3 11 0.4);

   prior r ~ uniform(2.6, 3.3);
   prior D ~ gamma(325.5, 0.0285714286);
   prior p ~ uniform(0.1, 0.4);

   initstate ID ~ degenerate(100);
   initstate IU ~ uniform(0.85*ID/p, 2*ID/p);
   initstate RC ~ degenerate(0);
   initstate S ~ degenerate(&N - ID - IU);
   initstate xi ~ poisson(r/(&N*D)*S*(ID+IU));

   state RC ~ degenerate(RC.l1 + 1/D*(ID.l1+IU.l1));
   state ID ~ degenerate((1-1/D)*ID.l1 + Y.l1);
   state IU ~ degenerate((1-1/D)*IU.l1 + xi.l1 - Y.l1);
   state S ~ degenerate(S.l1 - xi.l1);
   state xi ~ poisson(r/(&N*D)*S*(ID+IU));

   model y~normal(p*IU,ifn(IU>1,2*sqrt(p*(1-p)*IU),2*sqrt(p*(1-p))));
   learn nparticle=256 algorithm=sir method=PMCMC(nsample=25000
      nchain=8 nbi=5000 thin=5 seed=1
      outpost=mycas.out samset samsum
      sampler=rwm(ntu=250)
      statistics=all diagnostics=all);
run;

The summary tables of the model information are shown in Output 23.1.1.

The summary tables of the posterior samples are shown in Output 23.1.2.

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

The following statements load the smoothing estimates for each state variable into the Work library and then sort each data set by the time ID:

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

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

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

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

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

The following statements create a relevant data table, combine_S, in the Work library to record the smoothing estimates and actual values of the number of susceptible individuals. Then they draw a plot to compare the smoothing estimates to the actual values at each time step. The plot is shown in Output 23.1.4, where the shadow band indicates the range between the 2.5th and 97.5th percentiles of the smoothing estimates.

data combine_S;
   set smo_1 (keep=P_1 P_2 Q2);
   set sir_train (keep=t S);
run;

proc sgplot data=combine_S;
   band x=t upper=P_2 lower=P_1 /transparency=.3;
   series x=t y=S /legendlabel='Actual Value' lineattrs=(color=red);
   series x=t y=Q2 /legendlabel='Smoothing Median'
      lineattrs=(color=blue pattern=shortdash);
   yaxis label='S, Smoothing Median';
   xaxis label='Day';
run;

The following statements create a relevant data table, combine_ID, in the Work library to record the smoothing estimates and actual values of the number of documented active infectious individuals. Then they draw a plot to compare the smoothing estimates to the actual values at each time step. The plot is shown in Output 23.1.5, where the shadow band indicates the range between the 2.5th and 97.5th percentiles of the smoothing estimates.

data combine_ID;
   set smo_2 (keep=P_1 P_2 Q2);
   set sir_train (keep=t ID);
run;

proc sgplot data=combine_ID;
   band x=t upper=P_2 lower=P_1 /transparency=.3;
   series x=t y=ID /legendlabel='Actual Value' lineattrs=(color=red);
   series x=t y=Q2 /legendlabel='Smoothing Median'
      lineattrs=(color=blue pattern=shortdash);
   yaxis label='ID, Smoothing Median';
   xaxis label='Day';
run;

The following statements create a relevant data table, combine_IU, in the Work library to record the smoothing estimates and actual values of the number of undocumented active infectious individuals. Then they draw a plot to compare the smoothing estimates to the actual values at each time step. The plot is shown in Output 23.1.6, where the shadow band indicates the range between the 2.5th and 97.5th percentiles of the smoothing estimates.

data combine_IU;
   set smo_3 (keep=P_1 P_2 Q2);
   set sir_train (keep=t IU);
run;

proc sgplot data=combine_IU;
   band x=t upper=P_2 lower=P_1 /transparency=.3;
   series x=t y=IU /legendlabel='Actual Value' lineattrs=(color=red);
   series x=t y=Q2 /legendlabel='Smoothing Median'
      lineattrs=(color=blue pattern=shortdash);
   yaxis label='IU, Smoothing Median';
   xaxis label='Day';
run;

The following statements create a relevant data table, combine_xi, in the Work library to record the smoothing estimates and actual values of the number of new infectious individuals the next day. Then they draw a plot to compare the smoothing estimates to the actual values at each time step. The plot is shown in Output 23.1.7, where the shadow band indicates the range between the 2.5th and 97.5th percentiles of the smoothing estimates.

data combine_xi;
   set smo_4 (keep=P_1 P_2 Q2);
   set sir_train (keep=t xi);
run;

proc sgplot data=combine_xi;
   band x=t upper=P_2 lower=P_1 /transparency=.3;
   series x=t y=xi /legendlabel='Actual Value' lineattrs=(color=red);
   series x=t y=Q2 /legendlabel='Smoothing Median'
      lineattrs=(color=blue pattern=shortdash);
   yaxis label='xi, Smoothing';
   xaxis label='Day';
run;

The following statements create a relevant data table, combine_RC, in the Work library to record the smoothing estimates and actual values of the cumulative number of recovered individuals. Then they draw a plot to compare the smoothing estimates to the actual values at each time step. The plot is shown in Output 23.1.8, where the shadow band indicates the range between the 2.5th and 97.5th percentiles of the smoothing estimates.

data combine_RC;
   set smo_5 (keep=P_1 P_2 Q2);
   set sir_train (keep=t RC);
run;

proc sgplot data=combine_RC;
   band x=t upper=P_2 lower=P_1 /transparency=.3;
   series x=t y=RC /legendlabel='Actual Value' lineattrs=(color=red);
   series x=t y=Q2 /legendlabel='Smoothing Median'
      lineattrs=(color=blue pattern=shortdash);
   yaxis label='RC, Smoothing';
   xaxis label='Day';
run;

The median smoothing estimates of all the state variables are close to the actual values of the simulated data set. The following statements load the forecasting result for each state variable into the Work library and then sort each data set by the time ID:

data fc_4; set mycas.fc(where=(VariableID=4)); run;
proc sort data=fc_4; by _t; run;

data fc_6; set mycas.fc(where=(VariableID=6)); run;
proc sort data=fc_6; by _t; run;

The following statements create a relevant data table, combine_xi_forecast, in the Work library to record the forecasting estimates and actual values of the number of new infectious individuals the next day. Then they draw a plot to compare the forecasting estimates to the actual values at each time step. The plot is shown in Output 23.1.9, where the shadow band indicates the range between the 2.5th and 97.5th percentiles of the forecasting estimates.

data combine_xi_forecast;
   set fc_4 (keep=P_1 P_2 Q2);
   set sir_test (keep=t xi);
run;

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

The following statements create a relevant data table, combine_Y_forecast, in the Work library to record the forecasting estimates and actual values of the new confirmed cases. Then they draw a plot to compare the forecasting estimates to the actual values at each time step. The plot is shown in Output 23.1.10, where the shadow band indicates the range between the 2.5th and 97.5th percentiles of the forecasting estimates.

data combine_Y_forecast;
   set fc_6 (keep=P_1 P_2 Q2);
   set sir_test (keep=t Y);
run;

proc sgplot data=combine_Y_forecast;
   band x=t upper=P_2 lower=P_1 / transparency=.3;
   scatter x=t y=Y /legendlabel='Actual Value'
      markerattrs=(color=red size=6 symbol=circlefilled);
   series x=t y=Q2 /legendlabel='Forecasting Median'
      markers markerattrs=(color=blue size=6 symbol=triangle)
      lineattrs=(color=blue pattern=solid);
   yaxis label='Number of the New Observed Infections, Forecast';
   xaxis integer label='Day';
run;