#

# SSA
v2.1

#

#
Constructing the trajectory matrix, applying SVD decomposition, reconstruction
and forecasting

#

#
Dataset: Monthly US Imports from

#
Federal Reserve Bank of

#

#
Millions of dollars, not seasonally adjusted

#

#
Practical note: the exponential trend present in these data can be made linear
by a logarithmic transformation; note that,

# in
such a case, the analysis below should be change as the eigenvectors
controlling the trend and seasonal reconstruction

#
will be different.

#

#
Initialize

clear

ginit

clc

#
Read the data and divide by 1e+4

data_all
= read("data\usimpch9006.prn","double",198,1)/1e+4

#
Leave out the last 12 observations for forecasting

data
= data_all(1::186)

dataf
= data_all(187::end)

#
Construct the trajectory matrix and apply SVD decomposition

M =
80

T =
make_trajectory(data,M)

[L,Q]
= decompose_trajectory(T)

#
Examine the eigenvectors of the decomposition for identifying eigenvalues for
trend and seasonal reconstruction

#

gaddwin("Sample
Eigenvectors Corresponding to Trend")

gaddview(0.,0.5,1.,0.5)

gxaxis("linear",0,M,10)

gtitle("Eigenvectors
1, 4 and 5")

gcolor({"black","red","green"})

gplot(Q(:,{1,4,5}))

gaddview(0.,0.,1.,0.5)

gxaxis("linear",0,M,10)

gtitle("Eigenvectors
6, 9 and 15")

gcolor({"black","red","green"})

gplot(Q(:,{6,9,15}))

gaddwin("Sample
Eigenvectors Corresponding to Seasonality")

gaddview(0.,0.5,1.,0.5)

gxaxis("linear",0,M,10)

gtitle("Eigenvectors
2, 3 and 7")

gcolor({"black","red","green"})

gplot(Q(:,{2,3,7}))

gaddview(0.,0.,1.,0.5)

gxaxis("linear",0,M,10)

gtitle("Eigenvectors
8, 28 and 36")

gcolor({"black","red","green"})

gplot(Q(:,{8,28,36}))

# Now
in pairs

gaddwin("Eigenvector
pairs corresponding to seasonality")

gaddview(0.,0.,0.5,1.)

gtitle("Eigenvectors
2 and 3")

gcolor("red")

gplot(Q(:,2),Q(:,3))

gaddview(0.5,0.,0.5,1.)

gtitle("Eigenvectors
7 and 8")

gcolor("green")

gplot(Q(:,7),Q(:,8))

#
Reconstruct the trend using the eigenvalues corresponding to trend

hvec1
= { 1,4,5,6,9,15 }

[F1,f1]
= reconstruct_trajectory(T,L,Q,hvec1)

#
Plot the data and the fitted trend

gaddwin("Data
and Fitted Trend")

gxaxis("linear",0,190,10)

gxgrid("major")

gygrid("major")

gcolor({"black","red"})

gplot([data,f1])

# Now
reconstruct the seasonal component using the eigenvalues corresponding to
seasonal component

hvec2
= {2,3,7,8,10,11,12,13,14,15,16,17,18,19,20,21,22,23,28,34,35,36}

[F2,f2]
= reconstruct_trajectory(T,L,Q,hvec2)

#
Plot the seasonal component

gaddwin("Seasonal
Component")

gxaxis("linear",0,190,10)

gxgrid("major")

gygrid("major")

gcolor({"black","red"})

gplot(f2)

# If
the separation of the two components (trend and seasonal) is successful then
their weighted correlation should be small

wcorr
= weighted_correlation2(f1,f2,M)

print
"Weighted correlation between trend and seasonal components = ",wcorr

#
Note that we indeed achieve a very good separation...

# We
can now forecast the trend and seasonal components as well as the series itself
for the next 12 months

#

#
Forecast the trend

[g1,R,v2]
= forecast_trajectory(f1,Q(:,hvec1),12)

print
"Trend verticality coefficient = ",v2

gaddwin("Trend
Forecasts vs. Actual Values - Last 12 observations")

gtitle("Actuals
in black, forecasts in red")

gxaxis("linear",0,12,12)

gcolor({"black","red"})

gplot([dataf,g1(end-11::end)])

#
Forecast the seasonal

[g2,R,v2]
= forecast_trajectory(f2,Q(:,hvec2),12)

print
"Seasonal verticality coefficient = ",v2

gaddwin("Seasonal
Forecasts vs. Actual Values - Last 12 observations")

gtitle("Actuals
in black, forecasts in red")

gxaxis("linear",0,12,12)

gcolor({"black","red"})

gplot([dataf,g2(end-11::end)])

#
Now, add the forecasts - they track the actual values very well

g3 =
g1 + g2

gaddwin("Trend+Seasonal
Forecasts vs. Actual Values - Last 12 observations")

gtitle("Actuals
in black, forecasts in red")

gxaxis("linear",0,12,12)

gcolor({"black","red"})

gplot([dataf,g3(end-11::end)])

#
Compute the RMSE and MAE of the forecast errors

print
"RMSE = ",sqrt(colmse(dataf-g3(end-11::end)))

print
"MAE = ",colmad(dataf-g3(end-11::end))