计算差分方程的收敛点_时间序列分析第一章 差分方程(Simulation)

上一篇文章已经将差分方程的相关理论推导说明，本节将用代码将某些结论可视化。

一、一阶差分方程

1、递归法求解差分方程

给定y在t=-1期的值，以及w在t=0,1,2,...时期的值，则通过递归可以求出各时期的y值，代码如下

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
​
A = [0.5,1,2]  #设定Φ分别小于、等于、大于1
t = 1000  #给定总的时期为1000
​
for i,a in enumerate(A):
    y = [0]  #给定y在t=-1期的值
    print(f'动态模拟Φ={a}:')
    print('y-1=',y[-1])
    for j in range(0,t+1):
        w = np.random.randn()  #这里w是一个随机变量
        y_ = a*y[-1] + w      #由上一期计算当期
        y.append(y_)        #将当期加载y列表的最后
        print(f'y{j}=',y[-1])   #将动态模拟过程打印出来
    #这个条件判断是为了输出好看一点
    if i == 1:
        num = f'51{i+2}'
        #plt.subplot(int(num))
    elif i == 2:
        num = f'51{i+3}'
    else:
        num = f'51{i+1}'
    plt.subplot(int(num))
    plt.plot(range(-1,t+1),y,'-')   #将y的整个变化过程可视化
    plt.title(f'y(t) = {a}*y(t-1)+w(t)')

动态模拟过程结果如下

变化趋势可视化：

从图形中可以看到，Φ<1时，y围绕均值0波动，波动的幅度（方差）也基本没有改变，所以可以判断该时间序列是平稳的；Φ=1时，数据表现出明显的随机游走特征，有明显的符号集群现象，即前一期为‘+’，则后一期有很大概率也为‘+’；前一期为‘-’，则后一期很大概率也为‘-’，这是一个明显的单位根过程; Φ>1时，时间序列明显地发散。

2、动态乘子

在上一篇文章中，我们已经得出一阶差分方程的动态乘子为：

其大小只与时间间隔有关而与时期t无关。下面来看看对不同情况的Φ值，其动态乘子如何随时间间隔 j 变化。

%matplotlib inline
​
import numpy as np
import matplotlib.pyplot as plt
​
Φ = 0.9  #分别改变其值为-0.9,2，-2
n = 40
Φj = []
for j in range(1,n):
    Φj.append(Φ**j)
plt.bar(range(1,n),Φj)
plt.title(f'Φ={Φ}')

输出结果：

从图形可以看出，当|Φ|<1时，w对y的影响系数随时间间隔 j 的增大而逐渐减小，当 j 足够大时，动态乘子趋于零，0<Φ<1时按指数衰减，-1<Φ<0时振荡衰减；当|Φ|>1时，动态乘子随 j 的增大区域无穷。（注意为了图形效果明显我将n设为10）

3、长期效应

考虑w的某一冲击对y的累积效应：

%matplotlib inline
​
import numpy as np
import matplotlib.pyplot as plt
​
Φ = 0.9
n = 30
Φs = [0]
for j in range(0,n):
    Φs.append(Φs[-1]+Φ**j)
plt.bar(range(-1,n),Φs)
plt.title(f'Φ={Φ}')
plt.title(f'Φ={Φ}')

从图形可以看出，随着 j 的不断增大，累积效应最终会收敛于一个定值，理论告诉我们，这个定值为1/(1-Φ)=10。让我们取一个较大的n=100看看是否接近于10，从图形上可以看出基本为10。

二、P阶差分方程（二阶）

1、递归求解（动态模拟）

对于二阶差分方程，我们首先选取一些使得F的特征值的模分别大于、等于、小于1的特殊的Φ1和Φ2进行动态模拟，看看这些时间序列的图像趋势怎样。

使F的特征值为实数且模均小于1：Φ1=0.6，Φ2=0.2

使F的特征值为实数其中一个模大于1：Φ1=0.6，Φ2=3

使F的特征值为实数模均大于1：Φ1=0.6，Φ2=8

%matplotlib inline
​
import numpy as np
import matplotlib.pyplot as plt
​
Φ1 = 0.6   #分别将其值改为上面三种情况
Φ2 = 0.2
n = 500
y = [0,0]
for i in range(0,n):
    w = np.random.randn()
    y.append(Φ1*y[-1]+Φ2*y[-2]+w)
    print(f'y{i}=',y[-1])
plt.plot(range(0,n),y[2:],'-')
plt.title(f'y(t)={Φ1}*y(t-1)+{Φ2}*y(t-2)+w(t)')

输出结果

从图形可以看到，当Φ1，Φ2使得两个特征值均为绝对值小于1的实数时，二阶差分系统是稳定的，当Φ1，Φ2使得两个特征值之一为绝对值大于1的实数时，二阶差分方程不稳定（趋于无穷），并且λ越大，发散速度越快。

使F的特征值为共轭复数且模小于1：Φ1=0.2，Φ2=-0.2

使F的特征值为实数模均大于1：Φ1=0.5，Φ2=-1.5

可以看到当特征值的模小于1时，系统稳定，当特征值的模大于1时，系统是振荡发散的。

接下来我们直接用计算机模拟出稳定的Φ1及Φ2的取值区域。

%matplotlib inline
​
import numpy as np
import matplotlib.pyplot as plt
​
Φ1_s = np.linspace(-5,5,100)
Φ2_s = np.linspace(-5,5,100)
n = 100
#这里运用numpy的一些处理方式可能会使运行速度比三重循环快得多，但对numpy不熟
for Φ1 in Φ1_s:
    for Φ2 in Φ2_s:
        y = [0,0]
        s = 0
        for i in range(0,n):
            w = np.random.randn()
            y.append(Φ1*y[-1]+Φ2*y[-2]+w)
            if abs(y[-1])>100:    #在我们设定的初始值下，稳定系统的绝对值超过100的概率几乎为0
                s = 1
                break
        if s==0:
            plt.scatter(Φ1,Φ2)

输出结果

可以看到使二阶差分系统稳定的点分布的三角形和理论推导是非常吻合的，可以预想，如果代码中Φ1和Φ2取值足够密集（程序运行时间也足够长），那么这个三角形会被稳定点填满。

2、动态乘子

接下来我们看看不同Φ1和Φ2的动态乘子如何随 j 的增大而变化。

使F的特征值为实数且模均小于1：Φ1=0.6，Φ2=0.2

使F的特征值为实数其中一个模大于1：Φ1=0.6，Φ2=3

使F的特征值为实数模均大于1：Φ1=0.6，Φ2=8

%matplotlib inline
​
import numpy as np
import matplotlib.pyplot as plt
​
Φ1 = 0.6
Φ2 = 0.2
λ1 = (Φ1 + np.sqrt(Φ1**2+4*Φ2))/2
λ2 = (Φ1 - np.sqrt(Φ1**2+4*Φ2))/2
c1 = λ1/(λ1-λ2)
c2 = λ2/(λ2-λ1)
​
x =  range(0,30)
y = []
for j in x:
    y.append(c1*λ1**j+c2*λ2**j)
plt.bar(x,y)

结果

可以看到，两个特征值模均小于1时，动态乘子以指数速度衰减为0，而模大于1时，发散。

使F的特征值为共轭复数且模小于1：Φ1=0.5，Φ2=-0.8

使F的特征值为共轭复数且模大于1：Φ1=0.5，Φ2=-1.5

使F的特征值为共轭复数且模等于1：Φ1=0.5，Φ2=-1

%matplotlib inline
​
import numpy as np
import cmath
import matplotlib.pyplot as plt
​
Φ1 = 0.2
Φ2 = -0.2
λ1 = (Φ1 + cmath.sqrt(Φ1**2+4*Φ2))/2
λ2 = (Φ1 - cmath.sqrt(Φ1**2+4*Φ2))/2
c1 = λ1/(λ1-λ2)
c2 = λ2/(λ2-λ1)
​
x =  range(0,20)
y = []
for j in x:
    y.append(c1*λ1**j+c2*λ2**j)
plt.bar(x,y)
plt.title(f'Φ1={Φ1},Φ2={Φ2}')

模小于和大于1时，可以看到动态乘子分别为振荡衰减、发散，模为1时，动态乘子周期性的变化

3、长期效应

下面看看系统稳定时，w的冲击对y的累积效应。取Φ1=0.5，Φ2=-0.8

%matplotlib inline
​
import numpy as np
import cmath
import matplotlib.pyplot as plt
​
Φ1 = 0.5
Φ2 = -0.8
λ1 = (Φ1 + cmath.sqrt(Φ1**2+4*Φ2))/2
λ2 = (Φ1 - cmath.sqrt(Φ1**2+4*Φ2))/2
c1 = λ1/(λ1-λ2)
c2 = λ2/(λ2-λ1)
n = 50
​
x =  range(0,n)
y_s = [0]
y = []
for j in x:
    y.append(c1*λ1**j+c2*λ2**j)
    y_s.append(y_s[-1]+y[-1])
plt.bar(range(-1,n),y_s)
plt.title(f'Φ1={Φ1},Φ2={Φ2}')

结果

收敛值的理论结果为1/(1-Φ1-Φ2)=0.77，计算模拟结果也收敛于0.7左右。