Showing posts with label Programming. Show all posts
Showing posts with label Programming. Show all posts

Tuesday, April 08, 2025

Financial time series analysis with Matrix Pencil, the modern Prony's method - Python, Future prediction, Yahoo finance

 Originally Saturday, November 30, 2024

 The misprints of the elements in unitary-matrices in the recipe corrected on Tuesday 8th April 2025





The misprints  the elements in unitary-matrices in the recipe corrected on Tuesday 8th April 2025

 

 


 

Friday, March 28, 2025

Engle's ARCH motion prediction model with the simulation data with Python

 This is introduced in my cartoon video of Yukkuri Kaisetsu (Touhou Project fan-art):




Auto-Regressive Conditional Heteroscedasticity (ARCH)

ARCH was developed by an economist Robert F. Engle III having won the 2003 Nobel Memorial Prize in Economic Sciences for its achievement.

Dependent variable: the variance error terms of the first regression:
The explanatory variable X can be the lagged dependent variables and/or the other variables: 
Then, it find the coefficient γ of the lagged squared error terms with reference to the log-likelihood: 
 

Generalised Auto-Regressive Conditional Heteroscedasticity (GARCH)

GARCH assumes the variance of the error term symmetrically varies depending the average size of the error terms in pervious time steps. It adds the lagged variance on the explanatory variable of the second regression with reference to the log-likelihood for finding the coefficients γ and δ:

 

 Simulation Data with Python

The following exhibits display the simulation data evaluated using ARCH to illustrate how ARCH functions.

To facilitate the visual representation of this simulation, the most basic form of Engle's ARCH, as introduced in the Wikipedia entry below, has been implemented.

Ref: https://en.wikipedia.org/wiki/Autoregressive_conditional_heteroskedasticity

This simplified simulation demonstrates motion prediction for intercepting incoming flying projectiles with erratic movements, resembling the fluctuations of a stock price.

Following is my Python codes: 

Tuesday, November 26, 2024

Python, matplotlib.animation: Test showing changing values of Index Vs Stock

I am testing "matplotlib.animation" using the stock market index (NASDAQ in this case example) and the stock market price (Amazon in this case example)

# Parameter Adjustment

# Set Ticker Symbol
# e.g. Ticker="^GSPC" for S&P 500, Ticker="NDAQ" for NASDAQ, Ticker="TOPX" for TOPIX
Ticker_index="NDAQ" # Ticker sympol for the index
# Ticker="AMZN" for Amazon, Ticker="NTDOY" for NINTENDO
Ticker_stock="AMZN" # Ticker symbok for a stock

# Set the number of years to download data
No_Years_Data=1;

# Set the numbers of days for the future forecast
No_Days_FutureForecast=250*1


# For mathematics
import numpy as np
import math
import statistics

# For plotting
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import matplotlib.animation as animation
import itertools # for joining lists inside a list
from itertools import count

# Setting Plot Size
plt.rcParams["figure.figsize"] = (10,7)

# Calling stock values from Yahoo finance API
# Ref. https://www.kaggle.com/code/alessandrozanette/s-p500-analysis-using-yfinance-data
import yfinance as yf
%config InlineBackend.figure_format='retina'
import warnings
warnings.filterwarnings("ignore")

# Now, let's retrieve the data of the past x years using yfinance
Years="".join([str(No_Years_Data),'y']) # Set the number of years ! Adjust
Index = yf.Ticker(Ticker_index).history(period=Years)  # Ticker Symbokl ! Adjust !
Stock = yf.Ticker(Ticker_stock).history(period=Years)  # Ticker Symbokl ! Adjust !
# Let's take a look at the data
display(Index.tail())
display(Stock.tail())

# Change to list
# https://stackoverflow.com/questions/39597553/from-datetimeindex-to-list-of-times
Dates=Stock.index.tolist()
Dates=[Dates[t].strftime("%Y-%m-%d") for t in range(len(Dates))]
Values_Index=Index["Close"].tolist()
Values_Stock=Stock["Close"].tolist()

# Taking a difference of natural log of these values
y_1=np.diff(np.log(Index["Close"])).tolist()
y_2=np.diff(np.log(Stock["Close"])).tolist()



# Inline animations in Jupyter
# https://stackoverflow.com/questions/43445103/inline-animations-in-jupyter
epsilon=[] # storing errors
Window=30; # business days in 1.5 months
for b in range(math.floor(len(y_1)/Window)-1):
    st=0+Window*b; ed=st+Window; print(st,ed)

    plt.rcParams["animation.html"] = "jshtml"
    plt.rcParams['figure.dpi'] = 150

    fig, ax = plt.subplots()
    x_value = []
    y1_value = []
    y2_value = []
    epsilon_value=[]

    count_ = count();
    def animate(a):
        counts=next(count_)
        x_value=Dates[st+counts:ed+counts].copy()
        y1_value=y_1[st+counts:ed+counts].copy()
        y2_value=y_2[st+counts:ed+counts].copy()
        #epsilon_value=[abs((y1_value[w]-y2_value[w])/y2_value[w]) for w in range(len(y1_value))] # Relative
        epsilon_value=[abs((y1_value[w]-y2_value[w])) for w in range(len(y1_value))] # Absolute
        ax.cla()
        ax.plot(x_value,y1_value, label=Ticker_index, color='slategrey')
        ax.plot(x_value,y2_value, label=Ticker_stock, color='darkseagreen')
        plt.legend(loc="upper right")
        plt.xticks(rotation=90)
        plt.title(f'From {Dates[st]} to {Dates[ed+Window]}')
        ax.set_xlim(0,Window)
        epsilon.append(epsilon_value)
    Animation=animation.FuncAnimation(fig, animate, frames=Window, interval = 500)
    display(Animation)

# Separating figures
plt.show(block=False)

# Analysing errors
# joininig lists inside a list
epsilon=list(itertools.chain.from_iterable(epsilon))

print(statistics.mean(epsilon))
plt.figure(); plt.hist(epsilon); plt.show(block=False)



 

Monday, November 11, 2024

Normal Distribution Random Normal T-test Type I Error







I've edited and published Yukkuri Kaisetsu movie with my original stand picture of Statistical test based on Python programming.

 

Tuesday, December 12, 2023

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)


 

Monday, August 21, 2023

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