My codes of financial time series analysis and the AI art used in my Yukkuri Commentary movie posted on YouTube

This video presents our recent research on financial volatility forecasting, comparing classical econometric models, signal-processing techniques, and modern machine-learning approaches.

The study evaluates the forecasting performance, volatility prediction accuracy, computational efficiency, and interpretability of five different methods:

• GARCH
• SARIMA
• Matrix Pencil (MP)
• XGBoost
• Long Short-Term Memory (LSTM)

Using Bitcoin returns as a high-volatility financial asset, we investigate whether modern AI-based methods consistently outperform traditional forecasting techniques when accuracy, explainability, and computational cost are considered simultaneously.

Our results show that:

✓ SARIMA achieves the lowest return forecasting error.
✓ Matrix Pencil delivers comparable forecasting accuracy with extremely low computational cost.
✓ XGBoost provides a strong balance between accuracy and efficiency.
✓ LSTM exhibits relatively strong performance in tracking volatility patterns and market regime changes, despite higher forecasting errors and computational requirements.
✓ No single model dominates across all evaluation criteria.

The findings suggest that the Matrix Pencil method remains a competitive and highly interpretable alternative to black-box machine-learning models for financial forecasting applications. 

NicoNico: https://www.nicovideo.jp/watch/sm46580360 

Sakuya, Youmu, and Marisa challenge a unique cooking competition!
Using Bitcoin market data as ingredients, five dishes — GARCH, SARIMA, Matrix Pencil, XGBoost, and LSTM — to compete to predict future market movements ♬
This video explains the strengths and characteristics of traditional econometric models, signal-processing approaches, and modern AI techniques through a fun Yukkuri-style commentary.
Which model will create the best “forecasting recipe” for the future market?

Artwork including original character illustrations:
https://www.pixiv.net/en/tags/%E5%A6%...

Music, sound effects, and background materials used in this video are credited in the end credits.
#TouhouProject #Touhou #YukkuriCommentary #reimu #marisa #sakuya #youmu #remilia #yuyuko #Finance #Bitcoin #Cryptocurrency #QuantitativeFinance #FinancialEngineering #FinancialForecasting #VolatilityForecasting #Econometrics #ARIMA #SARIMA #GARCH #ArtificialIntelligence #MachineLearning #DeepLearning #DataScience #TimeSeriesAnalysis #ExplainableAI #XGBoost #LSTM #SignalProcessing #MatrixPencil #Prony #PronyMethod #Python #Research

This image is posted on my pixiv page: https://www.pixiv.net/en/artworks/147674029

My attempt is to generate an AI art with free of charge using ChatGPT. Many AI illustration tools charge fee and are not so flexible yet. I have experimented with generating an image based on my rough drawing together with the input script into ChatGPT.

This is an AI-synthetic illustration output by ChatGPT based on my rough drawing. The following is the input script. It is like asking the contemporary popular character designer to illustrate Touhou Project characters in his style. ☆彡

Redraw this illustration in a 1990s Japanese fantasy RPG illustration style with the scene shows three cooks working side by side in a cozy kitchen:

• on the left, a refined silver-haired maid in a blue-and-white uniform with a green ribbon, skillfully slicing vegetables with a kitchen knife

• in the center, a white-haired swordswoman in a green outfit, attentively stirring a pot of simmering food on the stove

• on the right, a blonde witch-like girl wearing an oversized black pointed hat and a black-and-white outfit, energetically stir-frying ingredients in a frying pan.

Please improve the line quality, proportions, and shading while preserving the original composition.

 

*  I have got to prepare for a cover image for a video I am currently creating to post YouTube. Furthermore, over here, it is too hot to be concentrated with my usual art works. .... seriously, it's tough enough to make me sleep quite a few hours! In an accommodation where I am currently living, there is no AC! My energy has recently been drained without noticing. Therefore, I just decided to spend my experiment in the AI output while my concentration and energy for my own artwork is sizably limited 

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

 

 

 

Lower Enveloping with the Different Interpolation Methods

Comparison 

Cost and Benefit Analysis - Basic simulation with SymPy - Econometrics with Python Programming

Cost and Benefit Analysis - Basic simulation with SymPy - Econometrics with Python Programming

 

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.

 

Modelling operational deposits balance and the prediction

 

Using MATLAB to create this moving window of the operational deposits balance based on Matrix Pencil (Prony-like method).

Leontief Model Simulation with Python Part 2: Experimental Attempt of Eigenvalue Problem when a matrix is singular i.e. non-invertible

 

This piece of my experimental attempt is based on the paper Singularity in the Discrete Dynamic Leontief Model by István Ábel1, Imre Dobos https://pp.bme.hu/so/article/view/8432/7719 .  This paper was accidentally found while searching for "Matrix Pencil for financial analysis".  Yet, this method is different from the Matrix Pencil (MP) method used in my current researching topic.  Nonetheless, this mathematical modelling applying the eigenvalue problem for the cross sectional econometric model is interesting enough to entice me. 

Honestly speaking, I have not fully understood their method behind because it is really off-topic from my current time series modelling. At the same time, I would like to keep this topic for my near future reference as I will be able to use my Linear Algebra application gained from my current research project. At least, I have attempted solving their mathematical modelling while learning by imitation. Therefore, it is likely to contain some error which I have to spend a sufficient amount of time for investigating and learning this realm of research. 

In terms of this example, the capital coefficient matrix C_i,i is singular containing a row containing zeros i.e. non-inversible. Therefore, the method of the eigenvalue problem is applied instead of using C_i,i^(-1).  Lambda (Greek letter) denotes the eigenvalue of the following equation. The error margin is kept below 1% or 5% at most.

My Python codes are displayed below:

 

# importing necessary tools

import matplotlib.pyplot as plt

import networkx as nx

import numpy as np

import pandas as pd

import random

import math

import statistics

import networkx as nx

import scipy.linalg

from scipy.stats import qmc

RefList=[

' https://en.wikipedia.org/wiki/Input%E2%80%93output_model'

,

' https://www.youtube.com/watch?v=KmVfmISjayA&t=134s'

,

' https://www.youtube.com/watch?v=z_9HwKet8G0&t=301s'

,

' https://www.sciencedirect.com/science/article/pii/S0895717710001093'

,

' https://pp.bme.hu/so/article/view/8432 https://pp.bme.hu/so/article/view/8432/7719'

]

print('Referring to the following articles and YouTube videos: \n')

for k in range(len(RefList)):

    print(RefList[k])

print(' \n')

print('x_i,t: the vector of output levels \n')

print('d_i,t: the vector of final demands (excluding investment) \n')

print('L: the Leontief input–output matrix \n')

print('C: the capital coefficient matrix \n')

print('C: The following equation is used when these matrices are non-singular i.e. inversible. \n')

# The following follows the example shown in the paper

# Index: Industrial Sectors

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

print(f'There are {len(Sctrs)} industrial sectors. \n')

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

Pdf_L = pd.DataFrame(data = L,index = Sctrs,columns = Sctrs)

display('The input–output matrix: L =',Pdf_L)

C = np.array([[0.3,0.4,0.45],[0,0,0],[0.6,0.8,0.9]])

Pdf_C = pd.DataFrame(data = C,index = Sctrs,columns = Sctrs)

display('The capital coefficient matrix: C =',Pdf_C)

# Identity Matrix

I=np.identity(len(C))

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

display(f'I: {len(C)} by {len(C[0])} identity matrix',Pdf_I)

print(' \n')

print('Leontief formula: x_i,t = L x_i,t + C _i,i (x_i,t+1 - x_i,t) \n')

print('Then, it is converted to:  x_i,t+1 = C_i,i^-1 {(I_i - L_i,i + C_i,i) x_i,t}    \n')

print('On the other hand, the capital coefficient matrix C_i,i is singular containing a row containing zeros i.e. non-inversible. \n ')

print('Therefore, the method of the eigenvalue problem is applied instead of using C_i,i^(-1) as follows. \n ')

print(f"Lambda (Greek letter) denotes the eigenvalue of the following equation (Ref. {RefList[4]} : ")

print('1. x_i,t+1 C_i,i = (I_i - L_i,i + C_i,i) x_i,t } \n')

print('2. Lambda x_i,t C_i,i = (I_i - L_i,i + C_i,i) x_i,t } \n')

print('3. Lambda C_i,i = (I_i - L_i,i + C_i,i) } \n')

# Eigenvalue finding

E= scipy.linalg.eig(C, (I - L + C))

for k in range(len(E[0])):

    print(round(E[0][k].real,18))

# Using the positive real number eigenvalue

EigVal=E[0].real

Where=(np.where(EigVal==max(EigVal)))[0][0]

Lambda=EigVal[Where]

print(f"The maximum eigenvalue used is Lambda= {Lambda} .")

# Corresponding Engenvector denoted as Rho (Greek letter)

ind=np.where(E[0] == Lambda)

ind[0][0]

Rho=E[1][:,ind[0][0]]

print(' \n')

# Let's denote the initial values of production of 3 sectors

x_0=np.array([392999.32,422999.27,446999.23])

Pdf_x_0 = pd.DataFrame(data = x_0,index = Sctrs)

display("Let's denote the initial values of production of 3 sectors: x_i,0 =",Pdf_x_0)

x_1=x_0/Lambda

Pdf_x_1 = pd.DataFrame(data = x_1,index = Sctrs)

display("Then, x_i,1 = x_i,0 / Lambda  =",Pdf_x_1)

Right=np.dot(L,x_0)+np.dot(C,x_1-x_0)

print(f"Then, the right hand side (I_i - L_i,i + C_i,i) x_i,t is {Right} , \n")

Left=x_0

print(f"and the left hand side (I_i - L_i,i + C_i,i) x_i,t is {Left}. \n")

print(f"The percentage difference between the left-hand side and the right-hand side is {(Left-Right)/Right}. \n")

 

ARCH, GARCH, and ARIMA

 

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 δ:

 

Auto-Regressive Integrated Moving-Average (ARIMA)

Another popular auto-regressive method (AR) + Integrated differences (I) + moving-averaging (MA). 

Ref.: https://people.duke.edu/~rnau/411arim.htm 

Integrated (I): Using the differentiating a non-stationary time series to attempt to transform it to a stationary series. 

Moving-Average (MR) model: Using the past forecasted error terms instead of the lags of the dependent variable. 

 Limitation of these AR methods

Ref. Modeling daily realized futures volatility with singular spectrum analysis    Dimitrios D. Thomakosa, Tao Wanga, Luc T. Willeb; ∗ Received 29 November 2001 https://doi.org/10.1016/S0378-4371(02)00845-2

Limited as the forecasting model: Whereas it provides the magnitude of the fluctuations e.g., volatility, their distribution is assumed to be symmetric. The error of the estimation expands for the future estimate

The coefficients of the lagged variable contain non-periodic noise terms, which disrupt the estimation of the important periodicities represented by these coefficients.

SSA and Prony-like methods can extract the cyclical physical components that contain important information about the periodicity and its magnitude by distinguishing them from the noise components.

SSA performs better than various auto-regressive methods and Hodrick-Prescott filter (Noise-filtering method) introduced in "Singular spectrum analysis for real-time financial cycles
measurement   Maximilien Coussin 1 https://doi.org/10.1016/j.jimonfin.2021.102532"

Python: Combination of Singular Spectrum Analysis (SSA) and Fast Fourier Transform (FFT)

 

 

This is my non-linear time series analysis combining Singular Spectrum Analysis (SSA) and Fast Fourier Transform (FFT). 

The dataset is the foreign exchange (FX) rate of Euro (EUR) based on US-Dollar (USD) obtained from Frankfurter App provided by the European Central Bank (ECB).

Implementing FFT for each reconstruction component (RC) of SSA gives us much more precise insights of the periodicity of the time series fluctuation. 

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

import requests

import json

# importing necessary tools

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 import stats

import statsmodels.formula.api as sm

import math

import cmath

from scipy.linalg import hankel

import scipy.odr

from scipy import odr

# For using Frankfurter app

import ast

# For graphing

import matplotlib.pyplot as plt

import matplotlib.gridspec as gridspec

import seaborn as sns

# Fast Fourier Transform (FFT) https://docs.scipy.org/doc/scipy/tutorial/fft.html

from scipy.fft import fft, ifft, fftfreq

# Showing the plots bigger

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

# defining imaginary number

i=complex(0,1)

i

## Downloading the time series data from Frankfurter App

# 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

#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 ky,vl in value.items():

        Rates_EURUSD.append(vl)

EURUSD=dict(zip(Dates_EURUSD,Rates_EURUSD))

print('First week:',list(EURUSD.keys())[0],',  Last week:',list(EURUSD.keys())[-1])

#Rates_EURUSD

### Graph output (1) ###

NW=6 # Number of windows to slice this time series dataset

LI=math.floor(len(EURUSD)/NW) # The length of interval

fig, axs = plt.subplots(NW)

for k in range(NW):

    axs[k].plot(Dates_EURUSD[k*LI:(k+1)*LI], Rates_EURUSD[k*LI:(k+1)*LI])

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

## Singular Spectrum Analysis (SSA)

print("The following refers to https://stats.stackexchange.com/questions/544105/singular-spectrum-analysis-and-their-eigentriplets")

print('This SSA method is almost identical to the one based on MATLAB "Singular Spectrum Analysis - Beginners guide" (MathWorks)')

print(' https://uk.mathworks.com/matlabcentral/fileexchange/58967-singular-spectrum-analysis-beginners-guide which I attempted to duplicate with Python. ')

# Re-defining the time series and its length

s=Rates_EURUSD;

N=len(Rates_EURUSD)

# Getting window length L and lagged length K

#L = 500-1

L = 100

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

       

### Graph output (2) ###

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

### Graph output (3) ###

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

# sample spacing

Days=N*7 # Number of days

Months=N*7//30 # Number of months

Years=N*7//360 # Number of years

T=1/Months # Defining the sample spacing step-size

yf = fft(s)

xf = fftfreq(N, T) #[:N//2]

# FFT for each RC

RC_fft=[]

RC_fft_power=[]

for l in range(L):

    RC_fft.append(fft(RC[:,l])) #[0:N//2])

    RC_fft_power.append(np.abs(RC_fft[l]))

# Denoise FFT of RC

dn_RC_fft=RC_fft.copy()

dn_RC_fft_power=RC_fft_power.copy()

for j in range(len(RC_fft)):

    for k in range(len(RC_fft[0])):

        DenoisingThreshhold=(sum(RC_fft_power[j])/len(RC_fft_power[j]))

        #DenoisingThreshhold=(max(RC_fft_power[j]))*0.90

        if RC_fft_power[j][k]<DenoisingThreshhold:

            dn_RC_fft_power[j][k]=0

            dn_RC_fft[j][k]=0

        else:

            continue

# Defining the cycle length

CycLen=(xf[0:N//2]) # Cycle length

#CycLen=CycLen[::-1]

CycLen

### Graph output (4) ###

# Showing some examples of the cycle lengths of the frequencies

fig, axs = plt.subplots(7)

for k in range(1,8):

    axs[k-1].stem(CycLen,dn_RC_fft_power[-k*10][N//2-1:-1]) #[0:N//2])

# Summing all the frequencies of RCs for the inverse FFT

tp_dn_RC_fft=list(map(list, zip(*dn_RC_fft)))

tp_dn_RC_fft_power=list(map(list, zip(*dn_RC_fft_power)))

sum_dn_RC_fft=[]

sum_dn_RC_fft_power=[]

for k in range(len(tp_dn_RC_fft)):

    sum_dn_RC_fft.append(sum(tp_dn_RC_fft[k]))

    sum_dn_RC_fft_power.append(sum(tp_dn_RC_fft_power[k]))

# ### Graph output (5) ### Uncomment the following to display

#plt.stem of the denoised power stemplot

plt.stem(CycLen, sum_dn_RC_fft_power[N//2-1:-1])

plt.grid()

plt.show()

### Graph output (6) ### Uncomment the following to display

# Inverse FFT to reconstruct the time series

plt.plot(s)

plt.plot(ifft(sum_dn_RC_fft).real)