Numerical scheme simulation of forward stochastic differential equations with Markov Switching

Preface

In the real world , Many phenomena or models have the characteristics of state switching , This phenomenon is described by Markov chain , It can well explain the random switching under different conditions , This leads to Markov switching of stochastic differential equations （SDEwMS） The rapid development in recent years . This stochastic mathematical model can not only be widely used in financial markets , And the same is true in the field of control engineering . However, the explicit solution of this kind of model is very difficult to find , In this article, we mainly establish the weak first order Euler and Milstein To solve such stochastic differential equations , Want to know the weak format of higher solution , You can refer to the following papers .

Yang Li, Taitao Feng, Yaolei Wang, and Yifei Xin. A high order accurate and effective scheme for solving Markovian switching stochastic models[J]. Mathematics, 2021, 9(6), 588.

One ： Stochastic differential equations with Markov switching （SDEwMS）

1： Markov chain knowledge

- Basic jump and Markov chain knowledge

Prior to “ Random knowledge ” In related articles , We mentioned in detail the knowledge of Markov chains and how to simulate discrete Markov chains , See the following article for details , I won't introduce , But the following format simulation will use this part of knowledge .
Discrete Markov chain simulation process

• Poisson random measure

Combine basic knowledge , We define $N(de,dt)$ For in space $\epsilon \times [0,T]$ The strength on the is $\lambda (de)dt$

Poisson random measure of , among $\lambda$ For in $\epsilon$ Upper $Lebanese$ measure ,

Yes $d{r}_t={\int }_{\epsilon }l(r_{t-,},e)N(de,dt)$ , among $l(i,e)=j-i,e\in {\Delta }_{ij};l(i,e)=0$, other ,

among ${\Delta }_{ij}$ It's length. Follow $q_{ij}$ Relevant inter cell , The definition is as follows ：

${\Delta }_{12}=[0,q_{12}),{\Delta }_{13}=[q_{12},q_{12}+q_{13}),...,{\Delta }_{1M}$

${\Delta }_{12}=[0,q_{12}),{\Delta }_{13}=[q_{12},q_{12}+q_{13}),...,$

${\Delta }_{1M}=\left[{\sum }_{j=2}^{M-1}{q}_{1j},{\sum }_{j=2}^M{q}_{1j}\right),$

${\Delta }_{21}=\left[{\sum }_{j=2}^M{q}_{1j},{\sum }_{j=2}^M{q}_{1j}+q_{21}\right),$

${\Delta }_{23}=\left[{\sum }_{j=2}^M{q}_{1j}+q_{21},{\sum }_{j=2}^M{q}_{1j}+q_{21}+q_{23}\right),..,$

${\Delta }_{2M}=\left[{\sum }_{j=2}^M{q}_{1j}+{\sum }_{j=1,j\ne 2}^{M-1}{q}_{2j},{\sum }_{j=2}^M{q}_{1j}+{\sum }_{j=1,j\ne 2}^M{q}_{2j}\right),...$

$=\left[{\sum }_{j=2}^{M-1}{q}_{1j},{\sum }_{j=2}^M{q}_{1j}\right),$

${\Delta }_{21}=\left[{\sum }_{j=2}^M{q}_{1j},{\sum }_{j=2}^M{q}_{1j}+q_{21}\right),$

${\Delta }_{23}=\left[{\sum }_{j=2}^M{q}_{1j}+q_{21},{\sum }_{j=2}^M{q}_{1j}+q_{21}+q_{23}\right),..,$

${\Delta }_{2M}=\left[{\sum }_{j=2}^M{q}_{1j}+{\sum }_{j=1,j\ne 2}^{M-1}{q}_{2j},{\sum }_{j=2}^M{q}_{1j}+{\sum }_{j=1,j\ne 2}^M{q}_{2j}\right),...$

among $M$ Is the length of the Markov chain .

2：SDEwMS Introduction to the structure of

In a complete probability space $\left(\Omega ,F,P\right)$ , With information flow ${\left\{F_t\right\}}_{t\ge 0}$ , commonly $SDEwMS$ It has the following differential structure ：

$d{X}_t=f(t,X_t,r_t^h)dt+g(t,X_t,r_t^h)d{W}_t,t\in (0,T],X_0=x_0$.

among $T$ Is a limited terminal moment , $r_t^h\in S=\{1,2,...,M\}$ Right continuous $Borel$ Measurable Markov chains .

$W_t={\left(W_t^1,...,W_t^d\right)}^T$ yes $d$ A Brownian motion in which the elements are independent of each other ,

$X_t={\left(X_t^1,...,X_t^q\right)}^T$ yes $q$ Dimensional random process , Drift term $f(\cdot ,\cdot ,\cdot ,):[0,T]\times R^q\times S\to R^q$ and

Diffusion term $g(\cdot ,\cdot ,\cdot ,):[0,T]\times R^q\times S\to R^{q\times d}$ Are all measurable functions , $x_0$ Is the initial value of the equation .

3： Existence and uniqueness theorem of solution

For a fixed constant $L$ And any finite terminal moment , $T>0,x,y\in R^q,\forall t\in [0,T]$, Drift coefficient

$f(\cdot ,\cdot ,\cdot ,)$ , Diffusion coefficient $g(\cdot ,\cdot ,\cdot ,)$ Satisfy $Lipschitz$ Conditions ：

$\left|f(t,x,i)-f(t,y,i)\right|\vee \left|g(t,x,i)-g(t,y,i)\right|\le L\left|x-y\right|,\left|f(t,x,i)-f(s,x,i)\right|\le L\left|t-s\right|$.
And linear growth conditions ：

$\left|f(t,x,i)\vee g(s,x,i)\right|\le M(1+\left|x\right|)$ ,

among $M=L\vee max\{\left|f(t,0,i)\vee g(t,0,i)\right|\}$.

So the equation $SDEwMS$ The above conditions are met , When the initial conditions $x_0\in F_0$ when , The equation has a unique solution

$\{X_t,t\in [0,T]\}$ .

For specific proof, please refer to relevant papers and literatures ！

Two ： Relevant Ito formula

Ito formula is a necessary theoretical knowledge for deriving relevant numerical schemes , So here is a detailed introduction .

1： Ito formula

For the above $SDEwMS$ equation , We remember $x\in R^q,\phi (t,x,i)\in C^{2,1}\left([0,T]\times R^q\times S;R\right)$ , Define operators

$L^0\phi ,L^j\phi ,L_e^{-1}\phi$ from $R^q\times S\to R:L_{\phi }^0=\frac{\partial \phi }{\partial t}+{\sum }_{k=1}^q\frac{\partial \phi }{\partial x_k}{f}_k+\frac{1}{2}{\sum }_{m=1}^d{\sum }_{l,k=1}^q\frac{{\partial }^2\phi }{\partial x_k\partial x_l}{g}_{l,m}{g}_{k,m}$

$+{\sum }_{j\in S}{q}_{ij}\phi$,

$L^j\phi ={\sum }_{k=1}^q\frac{\partial \phi }{\partial x_k}{g}_{k,j},L_e^{-1}\phi =\phi \left(t,x,i+l(i,e)\right)-\phi (t,x,i)$,

among $i\in S,i\le j\le d$.

We get the following ITO integral formula ：

$\phi \left(t_{n+1},X_{t_{n+1}},r_{t_{n+1}}\right)=\phi (t_n,X_{t_n},r_{t_n})+{\int }_{t_n}^{t_{n+1}}{L}^0\phi (s,X_s,r_s)ds+$

${\sum }_{j=1}^d{\int }_{t_n}^{t_{n+1}}{L}^j\phi (s,X_s,r_s)d{W}_s^j+{\int }_{t_n}^{t_{n+1}}{\int }_{\epsilon }{L}_e^{-1}\phi \left(s,X_s,r_s\right)\tilde{N}(de,ds)$,

here $\tilde{N}(de,ds)=N(de,ds)-\lambda (de)ds$ Is a compensated Poisson measure .

2： Ito unfolds

Combine the Ito formula operator mentioned above , To any $\phi \left(t_{n+1},X_{t_{n+1}},r_{t_{n+1}}\right)\in C^{2,1}(R^q\times S;R)$ , stay t_n Use Ito formula to get the following integral form ：

$\phi \left(t_{n+1},X_{t_{n+1}},r_{t_{n+1}}\right)=\phi \left(t_n,X_{t_n},r_{t_n}\right)+{\int }_{t_n}^{t_{n+1}}{L}^0\phi \left(s,X_s,r_s\right)ds+$

${\sum }_{j=1}^d{\int }_{t_n}^{t_{n+1}}{L}^j\phi \left(s,X_s,r_s\right)d{W}_s^j$,

In order to simplify the narration, ITO unfolds , We make $d=q=1$ , In the same way ：

$\phi \left(t_{n+1},X_{t_{n+1}},r_{t_{n+1}}\right)=\phi \left(t_n,X_{t_n},r_{t_n}\right)+{\int }_{t_n}^{t_{n+1}}{L}^0\phi \left(s,X_s,r_s\right)ds+{\int }_{t_n}^{t_{n+1}}{L}^1\phi \left(s,X_s,r_s\right)d{W}_s$,

We have to $L^0\phi \left(s,X_s,r_s\right)and{L}^1\phi \left(s,X_s,r_s\right)$ Use Ito formula to get ：

$\phi \left(t_{n+1},X_{t_{n+1}},r_{t_{n+1}}\right)=\phi \left(t_n,X_{t_n},r_{t_n}\right)+$

${\int }_{t_n}^{t_{n+1}}\left[L^0\phi \left(s,X_s,r_s\right)+{\int }_{t_n}^s{L}^0{L}^0\phi \left(\tau ,X_{\tau },r_{\tau }\right)d\tau +{\int }_{t_n}^s{L}^1{L}^0\phi \left(\tau ,X_{\tau },r_{\tau }\right)d{W}_{\tau }+{\int }_{t_n}^s{L}_e^{-1}{L}^0\phi \left(\tau ,X_{\tau },r_{\tau }\right)d{\tilde{N}}_{\tau }\right]ds+$

${\int }_{t_n}^{t_{n+1}}\left[L^1\phi \left(s,X_s,r_s\right)+{\int }_{t_n}^s{L}^0{L}^1\phi \left(\tau ,X_{\tau },r_{\tau }\right)d\tau +{\int }_{t_n}^s{L}^1{L}^1\phi \left(\tau ,X_{\tau },r_{\tau }\right)d{W}_{\tau }+{\int }_{t_n}^s{L}_e^{-1}{L}^1\phi \left(\tau ,X_{\tau },r_{\tau }\right)d\tilde{N}\tau \right]d{W}_s$.

In getting the above equation , if $\phi \left(s,X_s,r_s\right)$ Smooth enough , Let's continue to expand the multidimensional integrals involved above , Higher order numerical schemes can be obtained .

3、 ... and ： The numerical format

1： Definition of weak convergence

Given time division $0=t_0<t_1<...<t_{N-1}<t_N=T$ , hypothesis $\Delta t=t_{n+1}-t_n,\Delta W_n=W_{t_{n+1}}-$

$W_{t_n}(n=0,1,2,...,N)$, And there is a smooth function $G\in C_p^{2(\beta +1)}$ And the normal number $C$ Satisfy the inequality

$\left|E[G(X_t)-G(X^N)]\right|\le C\left(1+E[|X_0{|}^b]\right)(\Delta t{)}^{\beta }$,

among $\beta ,b\in N$ , be $X^N$ In order \beta Weakly converges to $X_T$. Of course, if there is weak convergence, there is strong convergence , When we don't take the expected absolute value with the whole trajectory error , Instead, take the expected absolute error of each point , Is our strong convergence error requirement , Higher order strongly convergent schemes are more complex , The conditions are more stringent , The order can be 0.5 The multiple relation of .

2： Weakly first order convergent scheme

We remember $f_k^i=f(t,X^n,i),g_k^i=g(t,X^n,i)$ , So there is the following format :
Set the initial conditions $X_0$ It is known that , about $0\le n\le N-1$ , about $m$ The second part of Vito Ito process $k$ Weight , We have the following format ：

1： Euler format （ weak 1.0 rank , strong 0.5 rank ）： $X_k^{n+1,i}=X_k^n+f_k^i\Delta t+g_k^i\Delta W_n$
2： Milstein format （ Strong and weak are 1.0 rank ）： $X_k^{n+1,i}=X_k^n+f_k^i\Delta t+g_k^i\Delta W_n+\frac{1}{2}{L}^1{g}_k^i\left((\Delta W_n{)}^2-\Delta t\right)$
3： Higher order convergent scheme

The form of the higher-order scheme is complex , And it is difficult to simulate and prove , Want to know the weak format of higher solution , You can refer to the following papers .

Yang Li, Taitao Feng, Yaolei Wang, and Yifei Xin. A high order accurate and effective scheme for solving Markovian switching stochastic models[J]. Mathematics, 2021, 9(6), 588.

Four ： Experimental simulation

• Markov chain code

import matplotlib.pyplot as plt
import numpy as np
import time
plt.rcParams['font.sans-serif']=['SimHei']## Chinese code scrambling ！
plt.rcParams['axes.unicode_minus']=False# The abscissa minus sign indicates a problem ！

start_time = time.time()

def markov_m(days,delta):
    original_value = 1  ##  Set initial state 
    r = np.array([[-1,0.6,0.4], [10 ,-5, 5], [3, 1, -4]])  ### The transfer rate matrix is given arbitrarily according to the definition 
    p = r * delta
    for i in range(r.shape[1]):
        p[i, i] = 1 + p[i, i]

    # P = np.exp(delta * r) - 1
    # for i in range(r.shape[1]):
    # P[i, i] = 1 + P[i, i]
    # print(P) ###delta The smaller it is ,P Follow p The closer the 

    U = np.random.uniform(0,1,days)
    q = np.zeros((1,days))

    for i in range(days):
        u = U[i]
        v_temp = original_value
        q[:,i] = v_temp
        original_value = 1
        s = 0
        while original_value <= p.shape[1] and u > s:
            s = s + p[v_temp - 1, original_value - 1]  ### Probability value superposition 
            original_value = original_value + 1  ### The matrix column index becomes larger 
        original_value = original_value - 1  ## Because of the beginning original_value =1, So we need to reduce 1
    return p,q.tolist()[0]

• SDEwMS Format simulation code

def schemes(a_ls,b_ls,names):

    time_ls = [16,32,64,128,256]
    tran_matx_ls = []
    state_vect_ls = []
    for i in time_ls:
        mres = markov_m(i, 0.1)## Call Markov chain simulation function 
        tran_matx = mres[0]
        state_vect = mres[1]
        print('\033[1;31m The length of time is %s The time transition probability matrix is ：\033[0m'%i)
        print(tran_matx)
        tran_matx_ls.append(tran_matx)
        print('\033[1;34m The length of time is %s The state vector is ：\033[0m'%i)
        print(state_vect)
        state_vect_ls.append(state_vect)

    logtimels = []
    logerrorls = []
    CRls = []
    msg = ''' 1：exp function  2：sin function  3： Of the original function 2.5 Power  4： Primitive function  '''
    print(msg)
    msg_ls = ['exp function ','sin function ',' Of the original function 2.5 Power ',' Primitive function ']
    ## Define initialization contents 
    choose = int(input(' Please select g(x) function ：'))
    echoes = 500
    error_ls = []
    value_matx = np.zeros((echoes,len(time_ls)))

    Euler = 'y + a * dt * y + b * y * dwt_temp'
    Milstein = Euler + ' + 0.5 * b * b * y * (dwt_temp * dwt_temp - dt)'
    schemes_ls = [Euler,Milstein]

    x0 = 1.5  ##SDE Initial value of true solution 
    for s,n in zip(schemes_ls,names):
        for e in range(echoes):
            for i in range(len(time_ls)):
                y = x0  ## Numerical format initial value 
                N = time_ls[i]
                dt = 1 / N
                yt = []
                dwt = np.random.normal(0,1,N) * np.sqrt(dt)
                nowstate_vect = state_vect_ls[i]

                ### The process of solving the true solution 
                v_ls = []
                vtemp = 0
                for j in range(N):## Cycle length 
                    a = a_ls[int(nowstate_vect[j] - 1)]## Take the corresponding mu coefficient 
                    b = b_ls[int(nowstate_vect[j] - 1)]## Take the corresponding sigm coefficient 
                    vtemp = vtemp + (a - 0.5 * b * b) * dt + b * dwt[j]## For discrete variable coefficients , Use the definition of integral to solve its true solution 
                    v_ls.append(vtemp)

                xT = ''## Allocate memory 
                finv = x0 * np.exp(v_ls[-1])
                if choose == 1:
                    xT = np.exp(finv)
                if choose == 2:
                    xT = np.sin(finv)
                if choose == 3:
                    xT = np.power(finv,2.5)
                if choose == 4:
                    xT = finv

                ## Numerical scheme solution 
                for j in range(N):
                    a = a_ls[int(nowstate_vect[j]-1)]## Take the corresponding mu coefficient 
                    b = b_ls[int(nowstate_vect[j]-1)]## Take the corresponding sigm coefficient 
                    dwt_temp = dwt[j]
                    y = eval(s)
                    if choose == 1:
                        yt.append(np.exp(y))
                    if choose == 2:
                        yt.append(np.sin(y))
                    if choose == 3:
                        yt.append(np.power(y,2.5))
                    if choose == 4:
                        yt.append(y)
                error_temp = yt[N - 1] - xT ## The difference between the end value of the true solution and the numerical solution 
                error_ls.append(error_temp)
            value_matx[e,:] = error_ls
            error_ls = []
        error_fin = np.abs(sum(value_matx) / echoes)## Weak convergence finally finds the expectation 
        log_time = np.log(1 / np.array(time_ls))
        log_error_fin = np.log(error_fin)
        CR = np.polyfit(log_time,log_error_fin,1)## First order polynomial fitting 
        print('\033[1;36m{0:*^80}\033[0m'.format(' The result of the calculation is '))
        print(' Format %s The absolute value of the error is expected to be ：%s'%(n,np.round(error_fin,6)))
        print(' Format %s The weak convergence order of is ：%s'%(n,round(CR[0],6)))
        logerrorls.append(np.round(log_error_fin,6))
        logtimels.append(log_time)
        CRls.append(round(CR[0],6))
    return logtimels,logerrorls,CRls,msg_ls[choose - 1]

statemu_ls = [1,0.8,0.9]## Drift coefficient 
statesigm_ls = [-0.1,0.2,0.05]## Diffusion coefficient 
scheme_names = ['Euler', 'Milstein']
R = schemes(statemu_ls,statesigm_ls,scheme_names)

chart 1

chart 2

Figure above 1 Sum graph 2 Is the result of the simulation experiment , In line with the theoretical results ！

• Drawing code

marker_ls = ['-*','-o']
def make_figure():
    plt.figure(figsize=(8, 6))
    for i,j,l,m in zip(R[0],R[1],scheme_names,marker_ls):
        plt.plot(i,j,m,label=l)
    plt.plot(np.array([-4,-3]),np.array([-4,-3]) + 2.5,'--',label='slope=1')
    plt.xlabel('stepsize logN')
    plt.ylabel('log(|E(Y_T-X^N)|)')
    plt.title(R[3])
    plt.legend()
    plt.show()
make_figure()

chart 3

5、 ... and ： summary

SDEwMS yes SDE A generalization of , More trial , But the corresponding numerical solution theory will become more complex . In recent years , Many people have also studied backward stochastic differential equations with Markov and made many theoretical breakthroughs .

Stochastic differential equation is a very deep and wide research field , Whether forward or backward , Whether with switching, jumping or reflection , Are worth studying , Interested friends can read more works in this field , The content of this article is only for the purpose of attracting jade ！