The synthetic data imitating the business cycle is my original although the functions of generating the complex time series refers to my supervisor's work.
The set of Python codes implementing Singular Spectrum Analysis is directly copied from "Singular spectrum analysis and their "eigentriplets" (Stack Exchange). This SSA method is almost identical to the one based on MATLAB "Singular Spectrum Analysis - Beginners guide" (MathWorks) which I attempted to duplicate with Python.
This set of Python codes succeeds in the reconstruction component (RC) as shown in the graph.
# importing necessary tools
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import seaborn as sns
import networkx as nx
import numpy as np
import pandas as pd
import random
import statistics
import scipy.linalg
from numpy import linalg
from scipy import stats
import statsmodels.formula.api as sm
import math
import cmath
from scipy.linalg import hankel
# Showing the plots bigger
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (20,10)
## Creating a synthetic time series imitating the business cycle
# Defining the function for the synthetic complex exponential
def syn_exp(amp,freq,x):
t=len(amp)
y=[]
for omg in range(len(x)):
yi=0
for ap in range(len(amp)):
yi+=amp[ap]*cmath.exp(freq[ap]*x[omg])
y.append(yi)
return y
# Defining complex number
i=complex(0,1)
i
# Creating the synthetic data imitating the business cycle
phi = [0.03, 3*i, 10*i, -1+30*i, -1-20*i, 0.01+5*i];
alpha = [4, 3, 2, 2, 2, 1];
# sampling rate
S = 100
# number of exponentials in the signal
n = 600
# f_syn: samples used by Prony's methhod
_2n_1=(2*n-1)+1
omega = [2*math.pi/S*(_0__2n_1) for _0__2n_1 in range(_2n_1)]
# Defininig the synthetic time series
f_syn = syn_exp(alpha, phi, omega)
# Defining the noise
noise=[np.random.normal(0) for k in range(len(f_syn))]
noise=noise/linalg.norm(noise)
# eps: noise level, also try 10^(-0)
eps = 10^(-1);
#eps = 0
# f = f_syn added with random noise
f = f_syn + eps*noise;
f
# Re-defining the time series and its length
s=f;
N=len(f)
## The SSA starts
print("The following refers to https://stats.stackexchange.com/questions/544105/singular-spectrum-analysis-and-their-eigentriplets")
print("The SSA method of this set of Python codes is identical to the embedded covariance method of https://uk.mathworks.com/matlabcentral/fileexchange/58967-singular-spectrum-analysis-beginners-guide")
# Getting window length L and lagged length K
L = 500-1
K = N - L + 1
# Constructing the time lagged Hankel matrix
X=np.zeros((K,L))
for m in range (0, L):
X[:,m] = s[m:K+m]
# Trajectory matrix
Cemb = np.dot(X.T,X)/K
# Eigen decomposition
eigenValues, eigenVectors = np.linalg.eig(Cemb)
idx = eigenValues.argsort()[::-1]
eigenValues = eigenValues[idx]
eigenVectors = eigenVectors[:,idx]
# Vectors of Principal Components
PC = np.dot(X,eigenVectors)
# Pre-allocating Reconstructed Component Matrix
RC = np.zeros((N, L))
# Reconstruct the elementary matrices without storing them
for k in range(L):
myBuf = np.outer(PC[:,k], eigenVectors[:,k].T)
myBuf = myBuf[::-1]
RC[:,k] = [myBuf.diagonal(j).mean()\
for j in range(-myBuf.shape[0]+1, myBuf.shape[1])]
# # First 6 RC
# fig, ax = plt.subplots(3,2)
# ax = ax.flatten()
# for k in range (0, 6):
# ax[k].plot(RC[:,k])
# ax[k].set_title(str(i))
# plt.tight_layout()
# Plotting a graph comparison
RawSSA=RC[:,0]
for k in range(L-1):
RawSSA=np.add(RawSSA,RC[:,k+1])
SCs=10 # Smoothed cycle components
SmoothedSSA=RC[:,0]
for k in range(SCs-1):
SmoothedSSA=np.add(SmoothedSSA,RC[:,k+1])
plt.subplot(2,1,1)
plt.title("Original vs. Reconstructed time series")
plt.plot(s[:])
plt.plot(RawSSA)
plt.subplot(2,1,2)
plt.title("Reconstructed raw vs. smoothed time series")
plt.plot(RawSSA)
plt.plot(SmoothedSSA)
