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 rho1 rho2/ lower=-0.9999 upper=0.9999;
parms msd1 msd2 esd1 esd2 / lower=1.e-6;
/* Specify the terms with lagged response variables */
deplag LagsForMink(LogMink) LogMink(lags=1) LogMusk(lags=1);
deplag LagsForMusk(LogMusk) LogMink(lags=1) LogMusk(lags=1);
/* Specify the bivariate trend */
array rwQ{2,2};
rwQ[1,1] = msd1*msd1; rwQ[1,2] = msd1*msd2*rho1;
rwQ[2,1] = rwQ[1,2]; rwQ[2,2]=msd2*msd2;
state alpha(2) type=RW W(I) cov(g)=(rwQ);
comp minkLevel = alpha[1];
comp muskLevel = alpha[2];
/* Specify the bivariate white noise */
array wnQ{2,2};
wnQ[1,1] = esd1*esd1; wnQ[1,2] = esd1*esd2*rho2;
wnQ[2,1] = wnQ[1,2]; wnQ[2,2]=esd2*esd2;
state error(2) type=WN cov(g)=(wnQ);
comp minkWn = error[1];
comp muskWn = error[2];
/* Specify the observation equation */
model LogMink = LagsForMink minkLevel minkWn;
model LogMusk = LagsForMusk muskLevel muskWn;
/* Specify an output data set to store component estimates */
output out=salesFor press;
run;
proc ssm data=furs;
/* Specify the ID variable */
id year interval=year;
/* Define parameters */
parms rho1/ lower=-0.9999 upper=0.9999;
parms msd1 msd2 esd1 esd2 / lower=1.e-6;
/* Specify the terms with lagged response variables */
deplag LagsForMink(LogMink) LogMusk(lags=1);
deplag LagsForMusk(LogMusk) LogMink(lags=1) LogMusk(lags=1);
/* Specify the bivariate trend */
array rwQ{2,2};
rwQ[1,1] = msd1*msd1; rwQ[1,2] = msd1*msd2*rho1;
rwQ[2,1] = rwQ[1,2]; rwQ[2,2]=msd2*msd2;
state alpha(2) type=RW cov(g)=(rwQ);
comp minkLevel = alpha[1];
comp muskLevel = alpha[2];
/* Specify the bivariate white noise */
array wnQ{2,2};
wnQ[1,1] = esd1*esd1; wnQ[1,2] = -esd1*esd2;
wnQ[2,1] = wnQ[1,2]; wnQ[2,2]=esd2*esd2;
state error(2) type=WN cov(g)=(wnQ);
comp minkWn = error[1];
comp muskWn = error[2];
/* Specify the observation equation */
model LogMink = LagsForMink minkLevel minkWn;
model LogMusk = LagsForMusk muskLevel muskWn;
/* Specify an output data set to store component estimates */
output out=salesFor press;
run;
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;