Novel econometric method: Short-Time Prony Transform with Python

 

The following data analysis is estimating the daily return (natural log difference) of the foreign exchange rate (FX) of Euro (EUR) based on US-Dollar (USD) from 1997 to 2023.

 

The data analysis introduced here is the brand-new approach to the financial economics which often involves non-linear variables with an exponentially expanding/decaying signal which the conventional econometric methods tend to struggle dealing with.

In terms of the non-linear signal processing like approach, many financial data analysts have adapted Fourier transfer. However, Fourier transfer has a critical disadvantage of dealing with the non-stable signal as shown in these graphs and failing to detect more than one signal frequencies close to each other.

This brand-new method adapts the signal processing method called Prony method invented by Baron Gaspard Clair François Marie Riche de Prony (22 July 1755 – 29 July 1839)

 

De Prony was the mathematician famous at the contemporary time period of when Fourier was alive. Fourier method has been popular for a long time period mainly because of the convenience with a limited computational capability. Nevertheless, the current advancement of computational mathematics has enabled to apply Prony method to the advanced mathematical data analyses.

Prony method is advantageous to detect various signal frequencies with their amplitude (coefficient) even for an unstable signal processing. Through the time length, it can describe various waving formations of a time series. For example, the seasonality and the event risk can be detected by observing these varying trends.

The conventional financial economic models are limited to finding one coefficient of the autocorrelation of either a variable or its error-terms/residuals. Furthermore, they are biased while dealing with a variable where no convergence can be proven. In contrast, this estimate based on Prony method figures out the entire dynamic movement of a variable. 

I have just started studying about this novel econometric method so that I admit that both my understanding and my implementation are still incomplete, and its improvement is quite likely to be required. At the same time, I am confident in having understood the core concept and the fundamentally required mathematical methods. Therefore, I have produced my Python codes implementing this data analysis method. 

Newly add function: Version with Short-Time Windows

The following Python codes decompose the time series of EUR/USD to transform into multiple time domain windows which break down data into years.

It allows to analyse the change in the signal frequencies and their amplitudes. It uses the Matrix Pencil Method (MPM) for finding the eigenvalues deriving the frequencies. 

It is more or less similar to Short-Time Fourier Transform (STFT). However, Fourier method can only depict stable signalling cycles. In contrast, this alternative model allows to describe the aperiodic cycle (expanding or contracting) of signals.

The linear algebra tool of solving the equation to find the amplitudes (alpha s) working in Jupyter Notebook somehow does not work in Google Colaboratory so that there are two sets of codes for each. 

 

# 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.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"] = (20,10)

# defining imaginary number

i=complex(0,1)

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

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

# Converting the rate to natural logarithm

ln_Rates_EURUSD=np.log(Rates_EURUSD)

ln_Rates_EURUSD

#plt.plot(ln_Rates_EURUSD)

### Defining the function for plotting a signal processing model ###

# Assuming the sampling rate used for reconstructing the signal processing model

M=100

# f_syn: samples used by Prony's methhod

def OMEGA(n):

    _2n_1=(2*n-1)+1

    omega = [2*math.pi/M*(_0__2n_1) for _0__2n_1 in range(_2n_1)]

    return omega

#omega=OMEGA(n)

#omega

# Defining the function generating a time series variable

# based on the given alpha s, phi s, and the time denoted as the small omega

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

# plot the real part of the signal

def plot_signal(amp,freq,x,StartDate):

    Mrend = 2*math.pi;

    Mint = math.pi/2;

    t=np.linspace(0, Mrend, 10000)

    M=len(t)

    x= [2*math.pi/M*(_0__2n_1) for _0__2n_1 in range(len(t))]

    F=syn_exp(amp,freq,x)

    F_real=[F[k].real for k in range(len(x))]

    plt.plot(F_real,label=StartDate)

    plt.title("Comparison of the signal frequencies of the time domain windows")

    plt.xlabel("'pi/2'_______________'pi'_______________'3pi/2'_______________'2pi'")

    plt.ylabel('Signal: real part')

    plt.legend()

    plt.show

    return F

# Plotting with defined alpha, phi, and omega (steps)

# plot_signal(alpha,phi,omega,0,-1,"Title")

###### Main Short-Time Prony Transform ######

# Defining the time series F

#F=ln_Rates_EURUSD

F=Rates_EURUSD

n=len(F)

N=len(F)

print("Length of data:",n)

# Number of the windows

Interval=52 # Dont' change!!

Cut=0

Windows=math.floor(n/Interval)

print('Number of the windows:',Windows)

DateStart_List=[]

DateEnd_List=[]

alpha_list=[]

phi_list=[]

for w in range(Windows-1):

#for w in range(1):

    #print('Date starts:',Dates_EURUSD[Cut],'  Date ends:',Dates_EURUSD[Cut+Interval])

    DateStart_List.append(Dates_EURUSD[Cut])

    DateEnd_List.append(Dates_EURUSD[Cut+Interval])

    f=F[Cut:Cut+Interval].copy()

    Cut=Cut+Interval+(1-math.ceil(((w+1)%5)/1000))

    # subtract the last week data once 5 years to make sure it starts from a beginning of the calendar

    #plt.plot(f)

    n=len(f) # number of exponentials

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

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

        x=H[h]

        #panda_df = pd.DataFrame(data = x)

        #display(panda_df)

    # Finding the Eigenvalue

    E= scipy.linalg.eig(H0, H1)[0]

    #panda_df_E = pd.DataFrame(data = E)

    #display(panda_df_E)

    V=np.transpose(np.vander(E, N=None, increasing=True))

    #panda_df_V = pd.DataFrame(data = V)

    #display(panda_df_V)

    # Calculating the coefficient alpha s with the matrix pencil method (MPM)

    # Compute the amplitudes (coefficients) by the mldivide V\f[d:]=A

    if len(f[d:])!=d:

        Adj=1

    else:

        Adj=0

    ######## Recovering alphas #######

    ###### For Jupyter Notebook ######

    #A=LA.solve(V.real, f[d+Adj:]) # Only pick the real part of V matrix! # Length?

    #alpha_rep=A

    #alpha=alpha_rep.tolist()

    ###### For google colaboratory ######

    A=LA.solve(V, f[d:])

    alpha_rep=A.real

    alpha_rep=A

    alpha=alpha_rep.tolist()

    # Recovering phi s by transforming the eigenvalues to their natural logarithm

    # phi is the power of the exponential which is the complex number

    # phi s represent the frequency of this signal

    # frequencies and damping factors

    phi_rep = [(cmath.log(E[t])) for t in range(len(E))] # .imag shows only the imaginary part

    phi=phi_rep

    # Removing relatively less significant alphas and phis

    # Adjust this number of terms to recover.

    NumTerms=9 # Define number of terms to recover # Adjust here!

    for l in range(len(alpha_rep)-NumTerms):

        abs_alpha=[abs(alpha[k]) for k in range(len(alpha))]

        #print(abs_alpha.index(min(abs_alpha)))

        #print(alpha[abs_alpha.index(min(abs_alpha))])

        alpha.pop(alpha.index(alpha[abs_alpha.index(min(abs_alpha))]))

        phi.pop(phi.index(phi[abs_alpha.index(min(abs_alpha))]))

        abs_alpha.pop(abs_alpha.index(min(abs_alpha)))

    #display(alpha)

    #display(phi)

    ## Showing the plots bigger

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

    ## Plotting with defined alpha, phi, and omega (steps)

    ##omega=OMEGA(n)

    ##print(n, len(omega))

    #plot_signal(alpha,phi,omega,Dates_EURUSD[Cut])

    # Add the list

    alpha_list.append(alpha)

    phi_list.append(phi)

###### Plotting each window of the Short-Time Prony Transform ######

YearInterval=4

NumInt=math.ceil(len(alpha_list)/YearInterval)

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

gs = gridspec.GridSpec(10, 10)

Cnt=0

for k in range(NumInt):

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

    for l in range(min(YearInterval,len(alpha_list)-YearInterval*k)):

        # Plotting with defined alpha, phi, and omega (steps)

        #omega=OMEGA(n)

        omega=OMEGA(N)

        plot_signal(alpha_list[Cnt],phi_list[Cnt],omega,DateStart_List[Cnt])

        Cnt=Cnt+1

 

Comparing the clusters with the different sampling rate following Shannon-Nyquist theorem

 

This Python programming tests how important to set the sampling rate denoted by M.

Shannon-Nyquist requirements claims that the sampling rate M has to be exact: either too small or too big. 

# importing necessary tools

import matplotlib.pyplot as plt

import matplotlib.gridspec as gridspec

import numpy as np

import pandas as pd

import random

import math

import cmath

# Defining the imaginary number

i=complex(0,1)

###### Drawing the full complex field circle with a different step size ######

# Definining the sampling rate

M=10

# Generating the list of the imaginary numbers

ImagList=[]

for j in range(1,M+1):

    ImagList.append((j/M)*i)

# Generating the exponentials of the imaginary numbers in ImagList

ExpList=[cmath.exp(2*math.pi*ImagList[j]) for j in range(len(ImagList))]

# Separating the exponentials into the real part and the imaginary part

ExpListR=[ExpList[j].real for j in range(len(ExpList))]

ExpListI=[ExpList[j].imag for j in range(len(ExpList))]

# # # Uncomment the following lines to show the circle graph

# # Plotting graph for a complex field circle

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

# plt.plot(ExpListR,ExpListI,'ro')

# plt.xlabel('Real field')

# plt.ylabel('Imaginary field')

# plt.show

###### ###### Comparison of the sampling rates ###### ######

###### ###### to test Shannon-Nyquist theorem  ###### ######

FreqSamples=[3*i,4*i,13*i]

SampleRateMs=[10,100,1000]

ColorList=['b','c','m']

Freq=[]

for k in range(len(SampleRateMs)):

    TempFreq=[FreqSamples[j]/SampleRateMs[k] for j in range(len(FreqSamples))]

    Freq.append(TempFreq)

ExpFreq=[]

for k in range(len(Freq)):

    TempExpFreq=[cmath.exp(2*math.pi*Freq[k][j]) for j in range(len(Freq[k]))]

    ExpFreq.append(TempExpFreq)

ExpFreqR=[]

ExpFreqI=[]

for k in range(len(ExpFreq)):

    BenchR=[1,0,-1,0]

    TempExpFreqR0=[ExpFreq[k][j].real for j in range(len(ExpFreq[k]))]

    TempExpFreqR1=BenchR+TempExpFreqR0

    ExpFreqR.append(TempExpFreqR1)

    BenchI=[0,1,0,-1]

    TempExpFreqI0=[ExpFreq[k][j].imag for j in range(len(ExpFreq[k]))]

    TempExpFreqI1=BenchI+TempExpFreqI0

    ExpFreqI.append(TempExpFreqI1)

panda_df = pd.DataFrame(data = ExpFreq)

display(panda_df)

# Compare plots with the different sampling rates

# Ref. https://stackoverflow.com/questions/37360568/python-organisation-of-3-subplots-with-matplotlib

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

gs = gridspec.GridSpec(3, 3)

for k in range(len(ExpFreq)):

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

    plt.plot(ExpFreqR[k],ExpFreqI[k],'ro', color=ColorList[k])

    plt.title(f'Sampling Rate M= {SampleRateMs[k]}')

    plt.xlabel('Real field')

    plt.ylabel('Imaginary field')  

Synthetic complex time series data: Python generator

The following Python code generates the synthetic time series based on real numbers and complex numbers.

This is useful to simulate the business cycles, the investment returns, and the other wave formed time series models of finance and macroeconomics.

You will need to adjust some parameters depending on your usage. 

Running this Python programme saves a csv file to Download folder in your PC but you may change owing to your preference.

# importing necessary tools

import os

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

import matplotlib.pyplot as plt

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

# Initializing real numbers

x = -1.0

y = 0.001

 

# converting x and y into complex number

z = complex(x,y);

# defining complex number

i=complex(0,1)

# printing phase of a complex number using phase()

print ("The phase of complex number is : ",end="")

print (cmath.phase(z))

# signal Business Cycle

phi = [0.1, 0.3*i, 1*i, -1+3*i, -1-2*i, 0.03+0.5*i, 11*i, -9*i, 30*i]

alpha = [4, 3, 2, 2, 2, 1, 1, 1, 0.1]

# Defining the function generating the complex time series

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

# sampling rate

M = 100

# number of exponentials in the signal

n = 600

# Generating the time trend for the complex time series

_2n_1=(2*n-1)+1

omega = [2*math.pi/M*(_0__2n_1) for _0__2n_1 in range(_2n_1)]

omega

# Generating the synthetic time series data of the phytical part (without noise)

f_syn = syn_exp(alpha, phi, omega)

# Generating the noise of the time trend

noise=[np.random.normal(0) for k in range(len(f_syn))]

noise=noise/LA.norm(noise)

# eps: noise level, also try 10^(-0)

eps =  10^(-1);

#eps = 0

# Combining the physical parts and the noise parts of the time series

# f = f_syn added with random noise

f = f_syn + eps*noise;

# plot the real part of the signal

def plot_signal(amp,freq,x):

    Mrend = 2*math.pi;

    Mint = math.pi/2;

    t=np.linspace(0, Mrend, 10000)

    M=len(t)

    x= [2*math.pi/M*(_0__2n_1) for _0__2n_1 in range(len(t))]

    F=syn_exp(amp,freq,x)

    F_real=[F[k].real for k in range(len(x))]

    plt.plot(F_real)

    plt.title('Signal: real part')

    plt.xlabel("'pi/2'__'pi'__'3pi/2'__'2pi'")

    plt.ylabel('')

    plt.show

plot_signal(alpha,phi,omega)

# Generating the time stamp

# ref: https://stackoverflow.com/questions/4479800/python-generate-dates-series

import datetime

stY=2023

stM=6

stD=1

fLen=len(f)

edY=stY+math.floor(fLen/365)+math.floor((stM+(fLen%365)/30)/12)

edM=math.floor((stM+(fLen%365)/30)%12)

edD=math.floor((stD+fLen%365%30))

dt = datetime.datetime(stY, stM, stD)

end = datetime.datetime(edY, edM, edD, 0, 0, 0)

step = datetime.timedelta(days=1)

TimeLine0 = []

while dt < end:

    TimeLine0.append(dt.strftime('%Y-%m-%d %H:%M:%S'))

    dt += step

TimeLine=TimeLine0[:fLen].copy()

len(TimeLine)

SensorId=[1]*fLen

len(f)

# Combining the time stamp and the time series data in Panda Dataframe

dict = {'Timestamp': TimeLine, 'SensorId': SensorId, 'Value': f, }

df = pd.DataFrame(dict)

df.tail()

# Download this time series data to your PC

FileLocation_Folder="C:\\Users\\OWNER\\Downloads" ## Change here for your PC!!!

Direct = (r"","\\ComplexTimeSeries.csv")

os.makedirs(FileLocation_Folder, exist_ok=True)  

df.to_csv(FileLocation_Folder.join(Direct))

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