/* This program demonstrates the classical decomposition of a time series: Y = T + C + S + I where Y = time series to be described T = Trend C = Cycle S = Seasonal I = Irregular */ /* This segment of the program generates a Deterministic Trend. */ data mc1; do t = 1 to 100; tr = 100.0 + 4.0*t; output; end; keep tr t; proc gplot data=mc1; symbol v=dot c=black i=join h=.8; title1 'Deterministic Trend Data'; title2 'X=Time Y=Trend Series'; axis1 order=(0 to 100 by 10) label=(f=duplex 'Time'); axis2 order=(-100 to 700 by 100) label=(f=duplex 'Trend Series'); plot tr*t / haxis=axis1 vaxis=axis2; run; /* This segment of the program generates a Deterministic Cycle. */ data mc2; do t = 1 to 100; c = 50*cos(3.1416*t/10); output; end; keep c t; proc gplot data=mc2; symbol v=dot c=black i=join h=.8; title1 'Cycle with a=50, w = 2pi/20 (period = 20 months), theta = 0'; title2 'X=Time Y=Cyclical Series'; axis1 order=(0 to 100 by 10) label=(f=duplex 'Time'); axis2 order=(-150 to 150 by 30) label=(f=duplex 'Cycle'); plot c*t / haxis=axis1 vaxis=axis2; run; /* This segment of the program generates a Deterministic Seasonal. The Seasonal Effects add to one. */ data seasonal; input mon s; datalines; 1 -50 2 -25 3 25 4 -25 5 -50 6 50 7 75 8 50 9 5 10 -25 11 -50 12 20 ; proc gplot data=seasonal; symbol v=dot c=black i=join h=.8; title1 'Seasonal Effects by Month'; title2 'X=Time Y=Seasonal Effects'; axis1 order=(1 to 12 by 1) label=(f=duplex 'Time'); axis2 order=(-150 to 150 by 30) label=(f=duplex 'Seasonal Effects'); plot s*mon / haxis=axis1 vaxis=axis2; run; /* Setting up the seasonal effects and the seasonal dummy variables. */ data mc3; do t = 1 to 100; month = mod(t,12); if month = 0 then month = 12; if month = 1 then s = -50; if month = 2 then s = -25; if month = 3 then s = 25; if month = 4 then s = -25; if month = 5 then s = -50; if month = 6 then s = 50; if month = 7 then s = 75; if month = 8 then s = 50; if month = 9 then s = 5; if month = 10 then s = -25; if month = 11 then s = -50; if month = 12 then s = 20; if month = 2 then dum2 = 1; else dum2 = 0; if month = 3 then dum3 = 1; else dum3 = 0; if month = 4 then dum4 = 1; else dum4 = 0; if month = 5 then dum5 = 1; else dum5 = 0; if month = 6 then dum6 = 1; else dum6 = 0; if month = 7 then dum7 = 1; else dum7 = 0; if month = 8 then dum8 = 1; else dum8 = 0; if month = 9 then dum9 = 1; else dum9 = 0; if month = 10 then dum10 = 1; else dum10 = 0; if month = 11 then dum11 = 1; else dum11 = 0; if month = 12 then dum12 = 1; else dum12 = 0; output; end; run; title 'Seasonal Effects and Seasonal Dummy Variables'; proc print data = mc3; /* This segment of the program generates the Irregular Component. */ data mc4; do t = 1 to 100; i = 10*rannor(t); output; end; run; proc gplot data=mc4; symbol v=dot c=black i=join h=.8; title1 'Irregular Component'; title2 'X=Time Y=Irregular Component'; axis1 order=(0 to 100 by 10) label=(f=duplex 'Time'); axis2 order=(-150 to 150 by 30) label=(f=duplex 'Irregular Component'); plot i*t / haxis=axis1 vaxis=axis2; run; data combine; merge mc1 mc2 mc3 mc4; data combine; set combine; y1 = tr; y2 = tr + c; y3 = tr + c + s; y4 = tr + c + s + i; run; title 'y1 = tr, y2 = tr + c, y3 = tr + c + s, y4 = tr + c + s + i'; proc print data = combine; var y1 y2 y3 y4; proc gplot data=combine; symbol v=dot c=black i=join h=.8; title1 'Deterministic Trend + Cycle'; title2 'X=Time Y=Trend + Cycle'; axis1 order=(0 to 100 by 10) label=(f=duplex 'Time'); axis2 order=(-100 to 700 by 50) label=(f=duplex 'Trend + Cycle'); plot y2*t / haxis=axis1 vaxis=axis2; run; proc gplot data=combine; symbol v=dot c=black i=join h=.8; title1 'Deterministic Trend + Cycle + Seasonal'; title2 'X=Time Y=Trend + Cycle + Seasonal'; axis1 order=(0 to 100 by 10) label=(f=duplex 'Time'); axis2 order=(-100 to 700 by 50) label=(f=duplex 'Trend + Cycle + Seasonal'); plot y3*t / haxis=axis1 vaxis=axis2; run; proc gplot data=combine; symbol v=dot c=black i=join h=.8; title1 'Trend + Cycle + Seasonal + Irregular'; title2 'X=Time Y=Trend + Cycle + Seasonal + Irregular'; axis1 order=(0 to 100 by 10) label=(f=duplex 'Time'); axis2 order=(-100 to 700 by 50) label=(f=duplex 'Tr + Cyc + Seas + Irreg'); plot y4*t / haxis=axis1 vaxis=axis2; run; /* Here we fit the classical decomposition to the data using Proc NLIN, the SAS nonlinear least squares procedure. Notice, that to avoid the "dummy variable trap," we have dropped the seasonal dummy variable for one of the months, here January. So the intercept (beta0) is interpreted as being the "January" intercept for the trend and all the other months' intercepts are obtained by adding the respective month's estimated gamma coefficient to the January intercept. For example, (beta0 + gamma2) is the February intercept for the trend. To determine the estimates of each month's seasonal effects, you simply take the mean of all of the months' intercepts and subtract this mean from the respective month's intercept. The estimated intercepts by month are as follows: Jan = 51.74 Feb = 51.74 + 21.36 = 73.1 Mar = 51.74 + 75.21 = 126.95 Apr = 51.74 + 23.99 = 75.73 May = 51.74 + 2.21 = 53.95 Jun = 51.74 + 93.26 = 145.0 Jul = 51.74 + 120.3 = 172.04 Aug = 51.74 + 100.6 = 152.34 Sep = 51.74 + 56.16 = 107.9 Oct = 51.74 + 24.69 = 76.43 Nov = 51.74 - 0.58 = 51.16 Dec = 51.74 + 69.97 = 121.71 . The mean of all of the months' intercepts is (51.74 + 73.1 + ... + 121.71)/12 = 1208.05/12 = 100.67. Therefore, the estimated monthly seasonal effects are as follows: Jan = 51.74 - 100.67 = -48.93 (Actual = -50) Feb = 73.1 - 100.67 = -27.57 (Actual = -25) Mar = 126.95 - 100.67 = 26.28 (Actual = 25) Apr = 75.73 - 100.67 = -24.94 (Actual = -25) May = 53.95 - 100.67 = -46.72 (Actual = -50) Jun = 145.0 - 100.67 = 44.33 (Actual = 50) Jul = 172.04 - 100.67 = 71.37 (Actual = 75) Aug = 152.34 - 100.67 = 51.67 (Actual = 50) Sep = 107.9 - 100.67 = 7.23 (Actual = 5) Oct = 76.43 - 100.67 = -27.24 (Actual = -25) Nov = 51.16 - 100.67 = -49.51 (Actual = -50) Dec = 121.71 - 100.67 = 21.04 (Actual = 20) . */ title 'Nonlinear Least Squares Estimation of Classical Decomposition Model'; proc nlin data=combine maxiter=5000; parms beta0 = 50 beta1 = 4 a = 50 w = 0.3146 theta = 0 gamma2 = 75 gamma3 = 125 gamma4 = 75 gamma5 = 50 gamma6 = 150 gamma7 = 175 gamma8 = 150 gamma9 = 105 gamma10 = 75 gamma11 = 50 gamma12 = 120; model y4 = beta0 + beta1*t + a*cos(w*t + theta) + gamma2*dum2 + gamma3*dum3 + gamma4*dum4 + gamma5*dum5 + gamma6*dum6 + gamma7*dum7 + gamma8*dum8 + gamma9*dum9 + gamma10*dum10 + gamma11*dum11 + gamma12*dum12; der.beta0 = 1; der.beta1 = t; der.a = cos(w*t + theta); der.w = -a*t*sin(w*t + theta); der.theta = -a*sin(w*t + theta); der.gamma2 = dum2; der.gamma3 = dum3; der.gamma4 = dum4; der.gamma5 = dum5; der.gamma6 = dum6; der.gamma7 = dum7; der.gamma8 = dum8; der.gamma9 = dum9; der.gamma10 = dum10; der.gamma11 = dum11; der.gamma12 = dum12; output out=yfit p=predict r=resid; run; /* Plot Fitted vs. Actual Values */ proc gplot data=yfit; symbol v=dot c=black i=join h=.8; title1 'Fitted Values of Classical Decomposition Model'; title2 'X=Time Y=Trend + Cycle + Seasonal + Irregular + = Predicted Value'; axis1 order=(0 to 100 by 10) label=(f=duplex 'Time'); axis2 order=(-100 to 700 by 50) label=(f=duplex 'Tr + Cyc + Seas + Irreg'); plot y4*t predict*t / overlay haxis=axis1 vaxis=axis2; run;