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)

 

Marisa Kirisame and Master Spark

 

Special tribute to the No.1 girl in 2023! I have painted Marisa Kirisame with Master Spark! My favourite special move in Touhou Project! 

Tool used: Clip Studio Paint 

These illustrations are posted in my Pixiv page so please visit https://www.pixiv.net/en/users/32334195!    

 

 

Cirno-chan and Dai-chan, Touhou Project

 

Tool used: Clip Studio Paint 

These illustrations are posted in my Pixiv page so please visit https://www.pixiv.net/en/users/32334195!    

 

Kochiya Sanae, Touhou Project

 

 

Tool used: Clip Studio Paint 

These illustrations are posted in my Pixiv page so please visit https://www.pixiv.net/en/users/32334195!   

I have actually put an effort for this posing of Sanae-chan: have taken this croquis technique for this posing! 

 

Leontief Model Simulation with Python Part 1: Simulating "Dynamic Input-Output Model" by RCD Sir with Python

I have simulated "Dynamic Input-Output Model" by RCD Sir shown in the following YouTube with Python.

This is my first step to experiment simulating Leontief's Input-Output econometric model: The first one is the basic model utilising the inverse matrix to find the output level.

The output is derived from the final demand multiplied by the inverse matrix of an identity matrix minus the coefficient matrix and the capital matrix.

A list of my Python code based on the mathematical modelling is as follows:

# importing necessary tools

import matplotlib.pyplot as plt

import numpy as np

import pandas as pd

import scipy.linalg

print('Simulating the example https://www.youtube.com/watch?v=KmVfmISjayA&t=134s')

Sctrs=['Sec1','Sec2','Sec3']

# A: Input Coefficient Matrix

A = np.array([[0.2, 0.3, 0.2],[0.4,0.1,0.2],[0.1,0.3,0.2]])

Pdf_A = pd.DataFrame(data = A,index = Sctrs,columns = Sctrs)

display('A: Input Coefficient Matrix',Pdf_A)

# B: Capital Coefficient Matrix

B = np.array([[0.1,0.2,0.1],[0.2,0.1,0.2],[0.1,0.2,0.1]])

Pdf_B = pd.DataFrame(data = B,index = Sctrs,columns = Sctrs)

display('B: Capital Coefficient Matrix',Pdf_B)

# G: Diagonal matrix of sector growth rates

G = np.array([[0.2,0,0],[0,0.2,0],[0,0,0.1]])

Pdf_G = pd.DataFrame(data = G,index = Sctrs,columns = Sctrs)

display('G: Diagonal matrix of sector growth rates',Pdf_G)

# F: Final Demand Vector

d = np.array([400,300,500])

Pdf_d = pd.DataFrame(data = d,index = Sctrs)

display('d: Final Demand Vector',Pdf_d)

# Identity Matrix

I=np.identity(len(A))

Pdf_I = pd.DataFrame(data = I,index = Sctrs)

display(f'I: {len(A)} by {len(A)} identity matrix',Pdf_I)

BG=np.dot(B,G)

Pdf_BG = pd.DataFrame(data = BG,index = Sctrs,columns = Sctrs)

display('BG is',Pdf_BG)

print('When the output vector is x \n x = Ax + BGx + d \n and therefore \n (I - A - BG)x = d .')

I_A_BG=I-A-BG

Pdf_I_A_BG = pd.DataFrame(data = I_A_BG,index = Sctrs,columns = Sctrs)

display('(I-A-BG) is',Pdf_I_A_BG)

print('Then \n x = (I - A - BG)^-1 d .')

Inv_I_A_BG=np.linalg.inv(I_A_BG)

Pdf_Inv_I_A_BG = pd.DataFrame(data = Inv_I_A_BG,index = Sctrs,columns = Sctrs)

display('(I-A-BG)^-1 is',Pdf_Inv_I_A_BG)

x= np.dot(Inv_I_A_BG,d)

Pdf_x = pd.DataFrame(data = x,index = Sctrs)

display('x = (I - A - BG)^-1 d =',Pdf_x)

My drawing of Kinnoull Hill, Perth, Scotland, UK