#

# SSA v2.1

#

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

#

# Dataset: Monthly US Imports from China, Jan-1990 to June-2006, FRED database series IMPCH

# Federal Reserve Bank of Saint Louis online database http://research.stlouisfed.org/fred2/

#

# 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))