Update record
1. Originally posted on 23rd July 2024
2. Title change from "Estimates based on the Matrix Pencil (MP) method for financial dataset (Prony method for finance)" to "Matrix Pencil (MP) method for financial time series estimate (Prony method for finance)" on 14th August 2024
3. Amended "N=len(f)" to "N=len(samples)" and "y = f" to "y = samples" in the function "MP(samples,nterms)" on 14th August 2024
4. Added the mathematical statement sheet on 14th August 2024
5. The dimension setting of Hankel matrix in the function plotting singular values (SVs) modified on 21 August 2024
This analysis is continued from Matrix Pencil (MP) methods for financial dataset (Prony method for finance). Last time, it demonstrated the full reconstruction. However, I have realised that it is not supposed to include the terms for the noisy terms (non-structural components). Matrix Pencil (MP) is essentially designed to only model the physical/structured components by filtering the rest noisy non-structured components.
The mathematical statement is as follows:
In order to detect the number of the terms representing the physical terms, Singular Value Decomposition (SVD) of the rectangular Hankel matrix is used. The singular values of the structural components are significantly higher than the non-structural counterparts. Furthermore, one set of the components is represented by the relatively identical singular values. Therefore, when the group changes from one to the other, there is a significant change from the previous values.
The number of the terms (NTERMS) is determined with reference to this SVD process. These parameters retrieved by MP successfully reconstructs a smoothed time series of Euro-Dollar FX rate as shown in the line plot. In addition, limiting the NTERMS to only the physical terms enables to stably estimate future values. This example exhibits 156 data points of the future.
The following is my code implemented:
## Importing libraries and setting up graphs and complex number
# Execute this cell to have access to the API of Frankfurter App
import requests
import json
import ast
# 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
# Showing the plots bigger
plt.rcParams["figure.figsize"] = (20,10)
# defining imaginary number
i=complex(0,1)
## Downloading the dataset: Euro-Dollar weekly data from European Central Bank (ECB)
# Connecting to the API https://www.frankfurter.app/docs/ to get the EUR/USD rates
url = "http://api.frankfurter.app/1999-01-01.."
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
# # Show the dataset plot
# Xticks=list(EURUSD.keys())
# f=list(EURUSD.values()) # Define the time series f
# xi = list(range(len(Xticks)))
# # plot the index for the x-values
# plt.plot(xi, f)
# plt.xlabel('Weeks')
# plt.ylabel('EUR/USD')
# plt.xticks(xi, Xticks, rotation=90)
# plt.title('Weekly EUR/USD FX rate')
# plt.legend()
# plt.show()
# Define basic parameters
f=list(EURUSD.values()) # Define the time series f
f_org=f.copy() # Saving the original data
n=len(f)
# Prony method
## Defininig the functions for the exponential analysis and reconstruction
# Definining the function of Matrix Pencil (MP) method,
# one of Prony-like methods.
def MP(samples,nterms):
# Firstly creating a rectangular hankel matrix of the direct imput of time series
n_samples=len(samples)
ncols_value=1/3
ncols=int(n_samples*ncols_value)
nrows=int(n_samples+1-ncols)
if nterms>math.floor(n_samples*ncols_value):
if min(ncols,nrows)==ncols: #Increase the number of columns
d=nterms-ncols
ncols=ncols+d
nrows=nrows-d
elif min(ncols,nrows)==nrows: #Increase the number of rows
d=nterms-nrows
ncols=ncols-d
nrows=nrows+d
H=hankel(samples)[0:nrows,0:ncols]
# Secondly, extracting the unitary matrix (left eigenvectors) of the hankel matrix
SVD_H=LA.svd(H, full_matrices=True, compute_uv=True, hermitian=False)
U=SVD_H[0]
U0=U[0:len(U)-1,0:nterms]
U1=U[1:len(U),0:nterms]
# Generating (snall) omega, the periodic time-steps.
n=math.floor(len(samples)/2)
_2n_1=(2*n-1)+1
M=1 # Sampling rate:
#According to Shannon-Nyquist sampling theorem, the sampling rate must be twice of the maximum signal length.
#However, M=1 the financial data not based on a particular signalling cycle model.
#omega = [2*math.pi/M*(_0__2n_1) for _0__2n_1 in range(_2n_1)] # Signal processing time steps
omega = np.arange(_2n_1) # using a simple discrete time steps Updated 17 07 2024
# Obtaining (big) Omega containing the frequencies and the time step
# Generalised engenvalues are decomposed by eig(U0\U1)
# https://fr.mathworks.com/help/matlab/ref/eig.html
OMEGA= scipy.linalg.eig(np.dot(np.linalg.pinv(U0),(U1)))[0]
# frequencies and damping factors
# phi_rec = [cmath.log(OMEGA[k])*M/(2*math.pi) for k in range(len(OMEGA))] # if using the signal process time step
phi_rec = [cmath.log(OMEGA[k]) for k in range(len(OMEGA))] # if using a simple discrete time steps Updated 17 07 2024
# Retrieving alpha
V=np.transpose(np.vander(OMEGA, N=len(samples), increasing=True))
# https://stackoverflow.com/questions/11479064/multiple-linear-regression-in-python
y = samples
X = V
coeffs = np.linalg.lstsq(X, y, rcond=None)[0]
alpha_rec=coeffs
alpha_rec=alpha_rec.tolist()
#shorting alphas and phs according to absolute values of alphas
abs_alpha=[abs(alpha_rec[k]) for k in range(len(alpha_rec))]
alpha_rec = [abs_alpha for _,abs_alpha in sorted(zip(abs_alpha,alpha_rec), reverse=True)]
phi_rec = [abs_alpha for _,abs_alpha in sorted(zip(abs_alpha,phi_rec), reverse=True)]
return alpha_rec,phi_rec
# Function for generating time series based on alphas and phis
def Exp(amp,freq,samples):
n=math.floor(len(samples)/2)
y=[]
# Generating x, the periodic time-steps.
_2n_1=(2*n-1)+1
omega = np.arange(_2n_1) # using a simple discrete time steps
for omg in range(len(omega)):
yi=0
for ap in range(len(amp)):
yi+=amp[ap]*cmath.exp(freq[ap]*omega[omg])
y.append(yi)
return y
# Displaying alphas and phis in a table
def PametersDisplay(AlphasPhis,nterms):
Header=['alphas','phis']
List=AlphasPhis
Table=dict(zip(Header,List))
panda_df = pd.DataFrame(data = Table)
pd.set_option('display.max_rows',nterms)
panda_df.head(nterms)
display(panda_df)
# Plotting two lines
def plot2(f,f_hat,Title,Xticks):
# generating the reconstructed data given alphas and phis
#print(Title)
xi = list(range(len(Xticks)))
plt.plot(f, label="original")
plt.plot(f_hat, label="reconstructed")
plt.legend(loc="upper right")
plt.title(Title)
plt.xlabel('Weeks')
plt.ylabel('EUR/USD')
plt.xticks(xi, Xticks, rotation=90)
plt.show
# Plotting two lines
def plot2_future(amp,freq,samples,NoWeeks,Title):
# Function for generating time series based on alphas and phis
N=len(samples)
N_new=N+NoWeeks
n=math.floor(N_new/2)
y=[]
# Generating x, the periodic time-steps.
_2n_1=(2*n-1)+1
omega = np.arange(_2n_1) # using a simple discrete time steps
for omg in range(len(omega)):
yi=0
for ap in range(len(amp)):
yi+=amp[ap]*cmath.exp(freq[ap]*omega[omg])
y.append(yi)
# generating the reconstructed data given alphas and phis
#print(Title)
plt.plot(samples, label="Original")
plt.plot(y, label="Reconstructed")
plt.legend(loc="upper right")
plt.title(Title)
plt.xlabel('Weeks')
plt.ylabel('EUR/USD')
plt.show
return y
# plot logged Singular Value Decomposition (SVD) and the log difference
# and returns the odd number indices of the terms representing the physical components
# for a reference to check how many exponential components can be grouped together.
# ref. https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.hankel.html
def plot_lnSVD_logdiff(samples):
# calculating logged SVs and their log difference
N=math.floor(len(samples))H=hankel(samples)[0:math.floor(N*2/3),0:math.floor(N*1/3)]
SVD=LA.svd(H, full_matrices=True, compute_uv=True, hermitian=False)[1]
ln_SVD=[math.log(SVD[k]) for k in range(len(SVD)) ]
d_ln_SVD=[abs(ln_SVD[k+1]-ln_SVD[k]) for k in range(len(SVD)-1) ]
# Plotting two line plots with each different scale
fig, ax1 = plt.subplots()
color = 'tab:blue'
ax1.set_xlabel('time (s)')
ax1.set_ylabel('exp', color=color)
ax1.plot(ln_SVD, color=color)
ax1.tick_params(axis='y', labelcolor=color)
ax2 = ax1.twinx() # instantiate a second Axes that shares the same x-axis
color = 'tab:red'
ax2.set_ylabel('sin', color=color) # we already handled the x-label with ax1
ax2.plot(d_ln_SVD, color=color)
ax2.tick_params(axis='y', labelcolor=color)
fig.tight_layout() # otherwise the right y-label is slightly clipped
plt.show()
# Selecting the indices of the significant terms representing the physical components
# The number of indices has to be odd as the first term with its index no 0 is the real term constant
# and the rest complex terms are the pairs
EnumSVD = [list(enumerate(d_ln_SVD))[k] for k in range(len(d_ln_SVD))]
TermIndices=[]
for k in range(math.floor(len(samples)/3*2/3)):
if EnumSVD[k][1]>statistics.mean(d_ln_SVD)*3: # !! ** Select the threshold of selecting terms ** !!
if EnumSVD[k][0]%2==0:
TermIndices.append(EnumSVD[k][0]+1)
else:
TermIndices.append(EnumSVD[k][0])
return list(set(TermIndices))
PhysicalTerms=plot_lnSVD_logdiff(f)
PhysicalTerms
## Set the number of terms (NTERMS) to reconstruct
NTERMS=PhysicalTerms[-1] # <-- Set up NTERMS by refering the SVD
# the number of the reconstruction terms i.e. the amplitude (alpha) + the frequency (phi)
# Adjust it until actually working for the reconstruction
# The fully fitting reconstruction requires the integer slightly smaller than a half of the time series length
# You may set up any integer for NTERMS e.g. small enough to extract only the trend
# Set the lag parameter e.g.4 to make the number of the reconstruction
## Output
# Showing the MP analysis and reconstruction result
print(f'NTERMS={NTERMS}')
Parameters=MP(f,NTERMS)
# generating the reconstructed data given alphas and phis
alpha_rec, phi_rec = Parameters
f_hat=Exp(alpha_rec,phi_rec,f)
f_hat=np.array(f_hat).real # Extracting only the real number
f_hat=f_hat.tolist() # Converting back to list from numpy
EDkeys=list(EURUSD.keys())
Tittle=(f'Estimate based on Euro vs. US-Dollar FX rate from {EDkeys[0]} to {EDkeys[-1]} Number of terms used {NTERMS}')
alpha_rec, phi_rec = Parameters
Future_Weeks=52*3 # ! ** Adjust **
f_future=plot2_future(alpha_rec, phi_rec,f,Future_Weeks,Tittle)
