Resources

Annual Sales of Mink and Muskrat Furs

/*--------------------------------------------------------------
                    SAS Sample Library

       Name: ssmex13.sas
       Description: Example program from SAS/ETS User's Guide,
              The SSM Procedure
       Title: Annual Sales of Mink and Muskrat Furs
       Product: SAS/ETS Software
        Keys: Multivariate Modeling, Mink and Muskrat Furs
        PROC: SSM
        Notes:

--------------------------------------------------------------*/

 title;

data furs;
   input Year  LogMink  LogMusk @@;
   year = mdy(1, 1, year);
   format year year.;
   datalines;
1850   10.2962   12.0752     1881   10.4957   13.6280
1851    9.9594   12.1791     1882   10.7277   13.8444
1852   10.1210   12.5863     1883   10.7687   13.8824
1853   10.1327   13.1102     1884   10.8646   13.8953
1854   10.6543   13.1466     1885   11.4167   13.6134
1855   10.8364   12.7531     1886   11.2451   12.7572
1856   11.0281   12.4638     1887   11.0714   12.8483
1857   11.0341   12.6191     1888   11.3269   12.7509
1858   11.2415   12.6556     1889   10.6152   12.3177
1859   11.0551   12.4461     1890   10.4800   12.6828
1860   10.7084   12.0855     1891   10.2914   13.2617
1861   10.3448   12.2357     1892   10.6517   13.6000
1862   10.8088   12.7230     1893   10.9711   13.7479
1863   10.6911   12.7857     1894   10.8359   13.3827
1864   11.0305   13.1417     1895   10.8452   13.4222
1865   11.0077   12.9441     1896   11.1595   13.6087
1866   10.8475   12.6786     1897   11.2433   13.2208
1867   10.9759   12.9292     1898   11.1620   13.2515
1868   11.2061   13.3344     1899   10.6416   13.4610
1869   11.2164   12.9096     1900   10.7359   13.5512
1870   10.2295   12.3556     1901   10.7751   13.7410
1871   10.3730   13.0036     1902   10.9616   14.3164
1872   10.5781   13.4657     1903   11.1057   14.2131
1873   10.7086   13.5514     1904   10.9091   13.7369
1874   11.0092   13.4180     1905   10.9330   13.8702
1875   11.1882   13.1689     1906   11.0030   13.4518
1876   11.2799   13.2765     1907   10.5756   12.9177
1877   11.2780   12.9880     1908    9.9774   12.5151
1878   11.3415   13.0940     1909    9.7902   12.6188
1879   11.0444   13.1218     1910    9.9891   13.5267
1880   10.4652   13.0775     1911   10.4045   13.7784
;
data add(keep=Year  LogMink  LogMusk);
   do i=1912 to 1923;
      year = mdy(1, 1, i);
      LogMink = .;
      LogMusk = .;
      output;
   end;
run;

 ODS NOPROCTITLE;
 proc append data=add base=furs; run;
     proc sort data=furs;
  by year;
 run;

proc ssm data=furs plots=residual;

   /* Specify the ID variable */
   id year interval=year;

   /* Define parameters */
   parms phi11 phi12 phi21 phi22;
   parms rho1 rho2/ lower=-0.9999 upper=0.9999;
   parms msd1 msd2 esd1 esd2 / lower=1.e-8;

   /* Define variables used to form T, W, and Q */
   array tmat{4,4};
   tmat[1,1] = 1; tmat[2,2] = 1;
   tmat[3,1] = 1; tmat[3,3] = phi11; tmat[3,4] = phi12;
   tmat[4,2] = 1; tmat[4,3] = phi21; tmat[4,4] = phi22;

   array qmat{4,4};
   qmat[1,1]=msd1*msd1; qmat[1,2]=msd1*msd2*rho1;
   qmat[1,3]=qmat[1,1]; qmat[1,4]=qmat[1,2];
   qmat[2,1]=qmat[1,2]; qmat[2,2]=msd2*msd2;
   qmat[2,3]=qmat[1,2]; qmat[2,4]=qmat[2,2];
   qmat[3,1]=qmat[1,3]; qmat[3,2]=qmat[2,3];
   qmat[3,3]=qmat[1,1]+esd1*esd1;
   qmat[3,4]=qmat[1,2]+esd1*esd2*rho2;
   qmat[4,1]=qmat[1,4]; qmat[4,2]=qmat[2,4];
   qmat[4,3]=qmat[3,4]; qmat[4,4]=qmat[2,2]+esd2*esd2;

   array wmat{4,2};
   wmat[1,1] = 1; wmat[1,2] = 0; wmat[2,1] = 0; wmat[2,2] = 1;
   wmat[3,1] = 1; wmat[3,2] = 0; wmat[4,1] = 0; wmat[4,2] = 1;

   /* Specify the state equation */
   state alpha(4) t(g)=(tmat) w(g)=(wmat) cov(g)=(qmat) a1(4) print=(T cov);

   /* Specify the components used in the observation equation  */
   comp minkVal = alpha[3];
   comp muskVal = alpha[4];

   /* Specify the observation equation */
   model LogMink = minkVal;
   model LogMusk = muskVal;

   /* Specify components for output purposes: elements of trend mu */
   comp minkLevel = alpha[1];
   comp muskLevel = alpha[2];

   /* Specify an output data set to store component estimates */
   output out=salesFor press;
run;

proc ssm data=furs plots=residual;

   /* Specify the ID variable */
   id year interval=year;

   /* Define parameters */
   parms phi12 phi21 phi22;
   parms rho1/ lower=-0.9999 upper=0.9999;
   parms msd1 msd2 esd1 esd2 / lower=1.e-8;
   phi11 = 0;

   /* Define variables used to form T, W and Q */
   array tmat{4,4};
   tmat[1,1] = 1;
   tmat[2,2] = 1;
   tmat[3,1] = 1; tmat[3,3] = phi11; tmat[3,4] = phi12;
   tmat[4,2] = 1; tmat[4,3] = phi21; tmat[4,4] = phi22;

   array qmat{4,4};
   qmat[1,1]=msd1*msd1; qmat[1,2]=msd1*msd2*rho1;
   qmat[1,3]=qmat[1,1]; qmat[1,4]=qmat[1,2];
   qmat[2,1]=qmat[1,2]; qmat[2,2]=msd2*msd2;
   qmat[2,3]=qmat[1,2]; qmat[2,4]=qmat[2,2];
   qmat[3,1]=qmat[1,3]; qmat[3,2]=qmat[2,3];
   qmat[3,3]=qmat[1,1]+esd1*esd1;
   qmat[3,4]=qmat[1,2]-esd1*esd2;
   qmat[4,1]=qmat[1,4]; qmat[4,2]=qmat[2,4];
   qmat[4,3]=qmat[3,4];
   qmat[4,4]=qmat[2,2]+esd2*esd2;

   /* Specify the state equation */
   state alpha(4) t(g)=(tmat) cov(g)=(qmat) a1(4) print=(T cov);

   /* Specify the components used in the observation equation  */
   comp minkVal = alpha[3];
   comp muskVal = alpha[4];

   /* Specify the observation equation */
   model LogMink = minkVal;
   model LogMusk = muskVal;

   /* Specify components for output purposes: elements of trend mu */
   comp minkLevel = alpha[1];
   comp muskLevel = alpha[2];
   comp totalLogSales = (0 0 1 1)*alpha;

   /* Specify an output data set to store component estimates */
   output out=salesFor press;
run;

proc sgplot data=salesFor;
  title "Forecasts for LogMusk and LogMink";
  title2 "(With Pointwise 95% Confidence Bands)";
  yaxis label="Mink and Muskrat Fur Sales (log scale)";
  band x=year lower=lower_LogMusk upper=upper_LogMusk ;
  scatter x=year y=LogMusk;
  series x=year y=forecast_LogMusk/ name="muskline"
     legendlabel="Forecasts for LogMusk";
  band x=year lower=lower_LogMink upper=upper_LogMink/
    fillattrs=GraphConfidence2 ;
  scatter x=year y=LogMink;
  series x=year y=forecast_LogMink / lineattrs=GraphFit2
     name="minkline" legendlabel="Forecasts for LogMink";

  refline '01jan1912'd / axis=x name="RefLine"
     lineattrs=GraphReference(pattern = Dash)
     legendlabel="Start of multistep forecasts";
  discreteLegend  "muskline" "minkline"  "RefLine" ;
run;
title;

proc sgplot data=salesFor;
    title "Smoothed Estimate of the Trend for LogMink";
    title2 "(With Pointwise 95% Confidence Bands)";
    band x=year lower=smoothed_lower_MinkLevel upper=smoothed_upper_MinkLevel;
    series x=year y=smoothed_MinkLevel;
    refline '01jan1912'd / axis=x name="RefLine"
        lineattrs=GraphReference(pattern = Dash)
        legendlabel="Start of multistep forecasts";
    DiscreteLegend  "RefLine" ;
 run;
 title;