SSA attempt functioning: The main SSA codes copied from StackExchange

  

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)

 

Attempt to use Singular Spectrum Analysis (SSA) for Frankfurter App with Python

Singular Spectrum Analysis (SSA) is nowadays a popular analysis method for financial time series. 

SSA is similar to Prony's method (the analogy would describe them as cousins) which can decompose a time series with non-periodic signalling cycles and a trend (Not a waving signal). 

Autoregressive method has a limitation of distinguishing the physical terms and the noise terms of the volatility, and it is not well suited for the future forecast for a non-linear time series.

In contrast,  SSA and Prony method can decompose such time series with a multi-dimensional decomposition process to offer more detailed estimates than the other conventional methods.

The part implementing SSA in the following Python codes refers to MATLAB codes displayed in https://uk.mathworks.com/matlabcentral/fileexchange/58967-singular-spectrum-analysis-beginners-guide.

This is an attempt to implement SSA for Frankfurter App with Python.

The following Python codes will be transferable to analysing the other financial data.

# Execute this cell to have access to the API of Frankfurter App

import requests

import json

# importing necessary tools

import matplotlib.pyplot as plt

import matplotlib.gridspec as gridspec

import networkx as nx

import numpy as np

import pandas as pd

import random

import statistics

import scipy.linalg

from numpy import linalg as LA

from scipy import stats

import statsmodels.formula.api as sm

import math

from scipy.linalg import hankel

# For plotting with multiple axis with different scales

# https://stackoverflow.com/questions/9103166/multiple-axis-in-matplotlib-with-different-scales

from mpl_toolkits.axes_grid1 import host_subplot

import mpl_toolkits.axisartist as AA

# Showing the plots bigger

plt.rcParams["figure.figsize"] = (18,9)

# Connecting to the API https://www.frankfurter.app/docs/ to get the EUR/USD rates

url = "http://api.frankfurter.app/1999-01-01.."

# Write your code here

resp = requests.get(url,{'to':'USD'})

#print(resp.url)

#print(resp.status_code)

#print("Type:", type(resp.content))

#print(resp.content)

print("Type:", type(resp))

# Preparing for the right encoding

import ast

byte_str = resp.content

dict_str = byte_str.decode("UTF-8")

mydata = ast.literal_eval(dict_str)

print("Type:", type(mydata))

#mydata #show the accessed data

EURUSD=(mydata["rates"])

#EURUSD

# Sorting the data into the Python dictionary form

Dates_EURUSD=[]

Rates_EURUSD=[]

for key, value in EURUSD.items():

    Dates_EURUSD.append(key)

    for k,v in value.items():

        Rates_EURUSD.append(v)

EURUSD=dict(zip(Dates_EURUSD,Rates_EURUSD))

#EURUSD

#Rates_EURUSD

print("Length of data:",len(Rates_EURUSD))

# Converting the rate to natural logarithm

ln_Rates_EURUSD=np.log(Rates_EURUSD)

ln_Rates_EURUSD

#plt.plot(ln_Rates_EURUSD)   # Uncomment to show plot!

# Creating the log-difference (Daily return) of EUR/USD

d_ln_Rates_EURUSD=np.diff(np.log(Rates_EURUSD))

len(d_ln_Rates_EURUSD)

#plt.plot(d_ln_Rates_EURUSD)   # Uncomment to show plot!

# Obtaining the length which is half of the total length of this variable.

# This is needed for creating the matrix of this variable i.e. n/2 * n/2 matrix

f=ln_Rates_EURUSD.copy()

#f=d_ln_Rates_EURUSD.copy()

# Creating Hankel matrix which is needed to

n = len(f)

d=(int(n/2-n%2*0.5)) # Subtracting the remainder if it is an odd number

d

H=hankel(f[0:-1])[0:d,0:d]

#print(f'\n\n \033[1m Hankel Matrix H \033[0m')

#panda_df = pd.DataFrame(data = H[h])

#display(panda_df)

# Finding the Generalised Eigenvalue of these Hankel matrices H0 and H1

print('\n\n Finding the Generalised Eigenvalue and Eigenvectors of these Hankel matrices H0 and H1 ')

Eig = scipy.linalg.eig(H)

Evalue= Eig[0]

Evectors = Eig[1:]

Evalue_real=Evalue.real

print('\n\n \033[1m Eigenvalue \033[0m')

panda_df = pd.DataFrame(data = Evalue)

display(panda_df)

print('\n\n \033[1m Eigenvectors \033[0m')

display(Eig[1:])

# Defnining the function plotting Singular Value Decomposition (SDV)

# This apprximate the number of the physical/true signals i.e. the number of true alpha s and phi s

def plot_SVD(A):

    SDV=LA.svd(A, full_matrices=True, compute_uv=True, hermitian=False)[1]

    Ln_SDV=[math.log(SDV[k]) for k in range(len(SDV)) ]

    plt.plot(Ln_SDV)

    plt.title('Singular Values')

    plt.xlabel('')

    plt.ylabel('Natural Log')

    plt.show

print('\n\n Singular Value decomposition (SVD) of the Hankel matrix H0')

# Refers to the plot to determine the number of terms (NTERMS) to recover

#plot_SVD(H) # Uncomment to show the SVD plot

# Determining NTERMS

SDV=LA.svd(H, full_matrices=True, compute_uv=True, hermitian=False)[1]

Ln_SDV=[math.log(SDV[k]) for k in range(len(SDV)) ]

NTERMS=0

SlopeSVD=1 # Adjust accoording to the data type and its log(SV) plot

for k in range(len(Ln_SDV)):

    if Ln_SDV[k]>SlopeSVD:

        NTERMS+=1

        #print(NTERMS)

        continue

    else:

        break

# Generating the reconstruction components for NTERMS

RCs = []

for k in range(NTERMS):

    RCtemp=Evectors[0][:,k]

    RCs.append(RCtemp)

# Showing each reconstruction components

print('\n\n \033[1m Reconstruction Components for NTERMS \033[0m')

panda_df = pd.DataFrame(data = RCs)

display(panda_df)

# For plotting with multiple axis with different scales

# https://stackoverflow.com/questions/9103166/multiple-axis-in-matplotlib-with-different-scales

plt.rcParams["figure.figsize"] = (7*NTERMS,30)

gs = gridspec.GridSpec(50,10)

for k in range(NTERMS):

    ax = plt.subplot(gs[k, 0]) # row k, col 0

    plt.plot(RCs[k])

# # Uncomment the entire following lines to show the graph of comparison

# # Summing all the reconstruction components

# Sum=0

# SRCs=[]

# for k in range(len(RCs[0])):

#     for l in range(len(RCs)):

#         Sum=Sum+RCs[l][k]

#     SRCs.append(Sum)

# len(SRCs)

# # Showing the sum of the reconstruction components

# plt.rcParams["figure.figsize"] = (20,10)

# plt.plot(SRCs)

# plt.xlabel("Time")

# plt.ylabel("Estimated price")

# plt.show

# # Comparing the original time series and the reconstructed time series

# # Cutting the length of f

# fcut=[(f[2*k-1]+f[2*k])/2 for k in range(math.floor(len(f)/2))]

# fcut=f

# # Streching the length of RCs

# StrSRCs=[0]*(len(SRCs)*2)

# for k in range(math.floor(len(f)/2)):

#     StrSRCs[2*k]=SRCs[k]

#     StrSRCs[2*k+1]=SRCs[k]

# len(StrSRCs)

# #StrSRCs=SRCs

# ## Plottting graph with multiple y-axis with different labels

# # https://stackoverflow.com/questions/9103166/multiple-axis-in-matplotlib-with-different-scales

# host = host_subplot(111, axes_class=AA.Axes)

# plt.subplots_adjust(right=0.75)

# par1 = host.twinx()

# par2 = host.twinx()

# YLabelL="f"

# YLabelR="rcs"

# VariableL=fcut

# VariableR=StrSRCs

# host.set_xlabel("Year")

# host.set_ylabel(YLabelL)

# par1.set_ylabel(YLabelR)

# p1, = host.plot( VariableL, label=YLabelL)

# p2, = par1.plot( VariableR, label=YLabelR)

# host.legend()

# host.axis["left"].label.set_color(p1.get_color())

# par1.axis["right"].label.set_color(p2.get_color())

# plt.draw()

# plt.xticks(rotation='vertical')

# plt.show()

a

Python: Hankel Matrices H0, H1, Generalised Eigenvalue of H0 and H1, and Singular Value Decomposition (SVD) used for the Matrix Pencil Method (MPM) and Singular Spectrum Analysis (SSA)

 

This process is the part of the Python code for Prony method introduced in the last post. The process of finding Hankel matrices, Generalised Eigenvalues, and Singular Value Decomposition (SVD) are the essential processes for both the modern Prony method (Matrix Pencil Method (MPM) and Singular Spectrum Analysis (SSA). 

 SSA is the non-parametric method which is said to be influenced by Prony method. For financial time series analyses, SSA is popularly used nowadays. SSA is better at capturing non-linearly time series with erratic fluctuations. SSA is popular due to the multi-resolution decomposition of time series into time trends and volatilities. 

The following Python codes operates the following actions:

- Downloads the foreign exchange rates of Euro based on USD

- Generating a time series variable EUR/USD

- Creating Hankel matrices H0 and H1 based on this time series

- Finding Generalised Eigenvalue and Eigenvectors of H0 and H1

- Singular Value Decomposition (SVD) of the Hankel matrix H0

 

# Execute this cell to have access to the API of Frankfurter App

import requests

import json

# importing necessary tools

import matplotlib.pyplot as plt

import networkx as nx

import numpy as np

import pandas as pd

import random

import statistics

import scipy.linalg

from numpy import linalg as LA

from scipy.stats import qmc

from scipy.optimize import newton

from scipy import stats

import statsmodels.formula.api as sm

import math

import cmath

from scipy.linalg import hankel

import scipy.odr

from scipy import odr

# Showing the plots bigger

plt.rcParams["figure.figsize"] = (12,6)

# defining complex number

i=complex(0,1)

i

# Connecting to the API https://www.frankfurter.app/docs/ to get the EUR/USD rates

url = "http://api.frankfurter.app/1999-01-01.."

# Write your code here

resp = requests.get(url,{'to':'USD'})

#print(resp.url)

#print(resp.status_code)

#print("Type:", type(resp.content))

#print(resp.content)

print("Type:", type(resp))

# Preparing for the right encoding

import ast

byte_str = resp.content

dict_str = byte_str.decode("UTF-8")

mydata = ast.literal_eval(dict_str)

print("Type:", type(mydata))

#mydata #show the accessed data

EURUSD=(mydata["rates"])

#EURUSD

# Sorting the data into the Python dictionary form

Dates_EURUSD=[]

Rates_EURUSD=[]

for key, value in EURUSD.items():

    Dates_EURUSD.append(key)

    for k,v in value.items():

        Rates_EURUSD.append(v)

EURUSD=dict(zip(Dates_EURUSD,Rates_EURUSD))

#EURUSD

#Rates_EURUSD

print("Length of data:",len(Rates_EURUSD))

# Converting the rate to natural logarithm

ln_Rates_EURUSD=np.log(Rates_EURUSD)

ln_Rates_EURUSD

#plt.plot(ln_Rates_EURUSD)   # Uncomment to show plot!

# Creating the log-difference (Daily return) of EUR/USD

d_ln_Rates_EURUSD=np.diff(np.log(Rates_EURUSD))

len(d_ln_Rates_EURUSD)

#plt.plot(d_ln_Rates_EURUSD)   # Uncomment to show plot!

# Obtaining the length which is half of the total length of this variable.

# This is needed for creating the matrix of this variable i.e. n/2 * n/2 matrix

X=d_ln_Rates_EURUSD.copy()

d=(int(len(X)/2-len(X)%2*0.5)) # Subtracting the remainder if it is an odd number

d

f=d_ln_Rates_EURUSD.copy()

n=len(f)

M=100  # sampling rate.

# M is supposed to be twice of the actual signal frequency.

# But, we have never known the true frequency of the financial data.

# So, it is assumed to be 100 because it seems to provide an appropriate answer

print(n) # Checking the total length of this variable.

# Creating Hankel matrix which is needed to

def mat_ge(samples):

    n_samples = len(samples)

    d = int(n_samples/2)

    H0=hankel(samples[0:-1])[0:d,0:d]

    H1=hankel(samples[1:])[0:d,0:d]

    return H0, H1

H0, H1 =mat_ge(f)

H=(H0,H1)

for h in range(len(H)):

    print(f'\n\n \033[1m Hankel Matrix H{h} \033[0m')

    panda_df = pd.DataFrame(data = H[h])

    display(panda_df)

# Finding the Generalised Eigenvalue of these Hankel matrices H0 and H1

print('\n\n Finding the Generalised Eigenvalue and Eigenvectors of these Hankel matrices H0 and H1 ')

Eig = scipy.linalg.eig(H0, H1)

Evalue= Eig[0]

Evectors = Eig[1:]

Evalue_real=Evalue.real

print('\n\n \033[1m Eigenvalue \033[0m')

panda_df = pd.DataFrame(data = Evalue)

display(panda_df)

print('\n\n \033[1m Eigenvectors \033[0m')

display(Eig[1:])

   

# Showing Verndermonde Matrix

Veig=np.transpose(np.vander(Evalue, N=None, increasing=True))

#print('\n\n Vandermonde Matrix of the eigenvalues')

#panda_df = pd.DataFrame(data = Veig)

#display(panda_df)

# Defnining the function plotting Singular Value Decomposition (SDV)

# This apprximate the number of the physical/true signals i.e. the number of true alpha s and phi s

def plot_SVD(A):

    SDV=LA.svd(A, full_matrices=True, compute_uv=True, hermitian=False)[1]

    Ln_SDV=[math.log(SDV[k]) for k in range(len(SDV)) ]

    plt.plot(Ln_SDV)

    plt.title('Singular Values')

    plt.xlabel('')

    plt.ylabel('Natural Log')

    plt.show

print('\n\n Singular Value decomposition (SVD) of the Hankel matrix H0')

plot_SVD(H0)

 

 

P.S. This following website shows the useful code to simulate Singular Spectrum Analysis (SSA) with Python

Introducing SSA for Time Series Decomposition
Python · MotionSense Dataset : Smartphone Sensor Data - HAR
https://www.kaggle.com/code/jdarcy/introducing-ssa-for-time-series-decomposition