はじめに

時系列解析とは
　時系列データを説明できる数理モデルを作ること

例えば東京の気温と雨量が以下のように得られたとして、

y_{気温} =3 \times x_{月}

のような、数式を立てるのが目的です。数式を得ることができれば

未来の予測
原理・本質の理解

をすることができます。

時系列モデルとは
前述の数理モデルには様々な方式がありますが、一般的に時系列データの持つ特徴は以下が挙げられます。

時系列データ = 短期の周期変動+周期変動+トレンド+外因性+ホワイトノイズ

短期の周期変動：前日の気温・数時間前の気温
周期変動：年単位の変動
トレンド：地球温暖化のような恒久的な増加
外因性：モデルに取り込めないような突発的なもの
ホワイトノイズ：温度計の計測誤差

これらの特徴を挙げられます。

短期の周期変動
前日の気温と今日の気温はそこまで差がないだろうと考えたとします。
前日の気温にノイズ$\epsilon$を足し他モデルを考えられます。

y_{気温} = y_{前日の気温} + \epsilon

こんな感じです。このモデルはARモデルと言われます。自己回帰モデル(Auto regression model)です。

全てノイズ
一人あたりの映画視聴時間を考えて見ましょう。
曜日の影響もなし、全てノイズである程度の時間を視聴しているみたいなモデルが立てられます

y_{視聴時間} = C + \epsilon

これがMAモデルです。

前日、面白かった映画を途中まで見ていたから今日も続きを見ようと考えた場合は

y_{視聴時間} = C + \epsilon +  \beta \epsilon_{前日}

みたいに立てられます。

季節性
気温は年単位で変動しています。月ごとにサンプルしている場合は12ヶ月前の気温データから、今月の気温が予想できます。これは一見ARモデルです。しかし仮に、次のようなモデルを考えた場合はどうなるでしょう。

今月の気温 = (去年の同月 \pm 1月　の気温の平均)

例えば、11月の気温を予測するために、去年の10~12月の気温を平均するモデルです。

これはARモデルでも表現できますが、12ヶ月分のジャンプ＋ARモデルと考えます。

同じ配信で別画面で開いて音声が少しズレている状態が季節性
状態の影響がジャンプして与えられるという事自体に意味があるときそれは季節性なのだ。

トレンド
地球温暖化により毎年0.1℃程度温度上昇しているとしましょう。これをモデル化すると

今年の気温 = 去年の気温 + 0.1 + ノイズ

になります。比例的に増加していきます。この性質をトレンドと呼びます。

これはARモデルですが、さきほどのARモデルと違い、時間が立つほど気温は増加し、最終的には爆発します。

ARモデルから導かれるAR特性方程式の解の絶対値が1以上であれば、値は爆発します。
この中で特に解の絶対値が1である統計モデルを単位根過程と呼びます。

ARIMAモデル
ここまでの統計モデルを合体したモデルは$ARIMA$モデルと呼ばれます。

モデルの推定
$ARIMA$モデルの同定の流れは

単位根検定（DF検定・ADF検定・PP検定）を行うことで差分過程が定常過程となるまで差分を繰り返す
差分回数は同じで複数の$ARIMA$統計モデルに対して実際の現象データに当てはめ、AIC情報量等を用いてARMA次数を決定
得られた定常過程とARMA次数に対して最尤法またはカルマンフィルタを使いパラメタ推定

モデルの検証

最後に時系列統計モデルの検証を行います。

今回のお話：DF検定

モデルの推定に必要になる単位根検定の基礎であるDF検定について説明します。

目的
以下の時系列モデルに従うとして、

x_{t} = \rho x_{t-1} + \epsilon_t

\epsilon_t \in N(0,\sigma^2 )

以下のような時系列データ$x_0 ,\cdots , x_n$に対して

$\rho = 1$か$\rho < 1$であるかを判別することがDF検定の目的です。

データ例

\rho = 1,\epsilon_t \in N(0,1)

\rho = 1,\epsilon_t \in N(0,10)

\rho = 0.9,\epsilon_t \in N(0,10)

\rho = 0.99,\epsilon_t \in N(0,10)

\rho = 0.99,\epsilon_t \in N(0,1)

\rho = 0.999,\epsilon_t \in N(0,1)

\rho = 0.9999,\epsilon_t \in N(0,1)

\rho = 0.9999,\epsilon_t \in N(0,10)

\rho = 1.001,\epsilon_t \in N(0,10)

\rho = 1.0001,\epsilon_t \in N(0,10)

わかること

$\rho > 1$だと数値が発散する。
$\rho < 1$だと分散は時刻によらなくなる
$\rho = 1$だと分散は時刻に比例する
$\rho = 1$に近づけば$\rho=1$のグラフに近くなる（当然か）

3種のグラフには見た目に明らかに違いがあったので人間であれば判定できそうです。

見た目で判断できるかどうか

上が

\rho = 0.99,\epsilon_t \in N(0,1)

下が

\rho = 1.00,\epsilon_t \in N(0,1)

になります。

$\rho$にほとんど差がない場合は見た目そこまで差がないです。

そこで、数量的に検証しようとするのがDF検定です！

1以上について
先に、前提として$\rho > 1$だと数値が発散するため、以後$\rho > 1$は考えません。

\rho = 1.001,\epsilon_t \in N(0,10)

本質的な違い

$\rho = 1$だと、分散は時刻に比例
$\rho < 1$だと、分散は一定

相関検定

$\rho < 1$を仮定する

x_{t} = \rho x_{t-1} + \epsilon_t

\epsilon_t \in N(0,\sigma^2 )

のモデルに対して、$\rho$は$x_t$と$x_{t-1}$の自己相関係数であると考えられます。

\rho = \frac{\sum x_{t-1} x_t}{\sum x_t^2 }

解法

x_{t} = \rho x_{t-1} + \epsilon_t

より

Var(x_{t} )  = \rho^2 Var(\rho x_{t-1})  + Var(\epsilon_t)

Var(x_{t} )  = \rho^2 Var(\rho x_{t-1})  + \sigma^2

\therefore Var(x_{t} )  = \frac{ \sigma^2 }{1- \rho^2}

値が収束する。

よって、分母はサンプル数をTとすると

\sum^T_{t=1} x_t^2 　= \frac{ \sigma^2 }{1- \rho^2} T

分子は

\sum x_{t-1} x_t = \sum x_{t-1} (\rho x_{t-1} + \epsilon_t)

 = \sum x_{t-1} \rho x_{t-1} + \sum x_{t-1} \epsilon_t

 = \rho \frac{ \sigma^2 }{1- \rho^2} T + \sum x_{t-1} \epsilon_t

ここまでは簡単にわかる。

無相関の確率変数の積の分散は、積になる

V[XY] = V[X] V[Y]

さらに、

\epsilon_t \in N(0,\sigma^2 )

Var(x_{t} )  = \frac{ \sigma^2 }{1- \rho^2}

がわかっているので。十分に収束していれば

Var(x_{t-1} )  = \frac{ \sigma^2 }{1- \rho^2}

と言える。よって

 x_{t-1}  \epsilon_t 　\sim  N ( 0, \frac{ \sigma^4 }{1- \rho^2} )

といえる。

 \sum x_{t-1} \epsilon_t 　\sim  N ( 0, \frac{ \sigma^4 }{1- \rho^2} T  )

\rho = \frac{\sum x_{t-1} x_t}{\sum x_t^2 }

にいままでの結果を代入すると

\frac{\sum x_{t-1} x_t}{\sum x_t^2 } = \frac{\rho \frac{ \sigma^2 }{1- \rho^2} T + N ( 0, \frac{ \sigma^4 }{1- \rho^2} T  )}{ \frac{ \sigma^2 }{1- \rho^2} T }=  \rho + N(0,(1-\rho^2)/T)

が得られる。

検定
帰無仮説として$\rho = 0$を立てれば、統計量

\frac{1}{ \sqrt{T}}\frac{\sum x_{t-1} x_t}{\sum x_t^2 }  \sim   N(0,1)

であることから$\rho = 0$検定ができる。ただし、分布の収束を前提としている箇所が複数出てきているので、サンプル数$T$の大きさによっては使えない。単位根過程$\rho=1$に至っては分散が発散し、収束しないため使えない。

回帰係数検定

次に、回帰係数を推定し$\hat{\rho}$を計算したときの、回帰係数$\hat{\rho}$の分布と検定を考えてみよう。

推定誤差$e_t = y_t - \hat{\rho} y_{t-1} $の分布はどうなるか。

\hat{\rho}
=
\frac{\sum x_{t-1} x_t}{\sum x_t^2 } = \frac{\rho \frac{ \sigma^2 }{1- \rho^2} T + N ( 0, \frac{ \sigma^4 }{1- \rho^2} T  )}{ \frac{ \sigma^2 }{1- \rho^2} T }=  \rho + N(0,(1-\rho^2)/T)

より

\rho - \hat{\rho} = N(0,(1-\rho^2)/T)

となる。

e_t = y_t  - \hat{\rho} y_{t-1}  =  \rho y_{t-1} + \epsilon_t   - \hat{\rho} y_{t-1}  = 
(\rho - \hat{\rho} ) y_{t-1}  + \epsilon_t

よって

e_t \sim N(0,\frac{(1-\rho^2) y_{t-1}^2}{T} + \sigma^2 )

とわかる。

 \sum e_t^2 　\rightarrow \sum \frac{(1-\rho^2) y_{t-1}^2}{T} + \sigma^2 T

に収束すると考えられる。新たに$\sigma_{\rho}^2$を定義する。

\sigma_{\rho}^2 =  \frac{\sum e_t^2}{\sum y_{t-1}^2 } 　\rightarrow  \frac{(1-\rho^2) }{T} + \frac{\sigma^2 T}{\sum y_{t-1}^2 } \rightarrow  \frac{(1-\rho^2) }{T}

第二項は$\sigma^2 T >> \sum y_t $と考えられれれば$0$に収束することを使った。

\rho - \hat{\rho} = N(0,(1-\rho^2)/T)

と比較しながら、統計量$t$を次のように定義する。

t = \frac{\rho - \hat{\rho}}{\sigma_{\rho}} \sim N(0,1)

都合のいいことに$N(0,1)$に収束するのである。

言いたいこととしては、

\sigma_{\rho}^2 =  \frac{\sum e_t^2}{\sum y_{t-1}^2 }

を計算すれば、$\hat{\rho}$の推定誤差分散になるということである。

まとめ

\hat{\rho} = \frac{\sum x_{t-1} x_t}{\sum x_t^2 }

\rho - \hat{\rho} \sim  N(0,\frac{\sum e_t^2}{\sum y_{t-1}^2 } )

DF検定

本題に入ろう。
ここまでの議論は分散の収束を原理としている。単位根過程において分散は発散するため、前述2つの方法は使えない（使えないならなんで説明したねんー＞前述の統計量を単位根過程に適用した場合を考えるため）。ここからは、単位根過程を仮定した場合の分布・検定を考えていく。

目的としては、帰無仮説として$\rho = 1$とした場合に使える統計量を考えることだ。

ブラウン運動

DF検定では単位根過程のうちブラウン運動のみを対象とした検定である。その他の単位根過程の検定はDF検定を拡張した、ADF検定を用いて行うことになる。今回はDF検定のみを考えよう。まず、ブラウン運動は次のように定義される。

x_0 = 0

x_{t_1} = x_{t_2} + N(0,\sigma^2 (t_1 - t_2))

移動量は経過時間に比例して増加するという考えだ。先程の$\rho = 1$と同じモデルである。

シミュレート時間間隔$dt$の大小で時系列データの見た目も変わる。

$x_t$のヒストグラムを作ると

こんな感じで正規分布に従う。シミュレーション上$dt$の単位が入ることになるが、どんな$dt$でもこの正規分布は変わらない。

シミュレーションは以下のようにできる。

t = 0
x = 0

while True:
    
    dx = random.gauss(0,sigma*sq_dt)
    x = x + dx    

    if t > T:
        break
    else:
        t = t + dt

単純にシミュレーション間隔$dt$ごとにホワイトノイズを足すだけである。

グラフ表示まで全体のプログラム以下のようになる。

プログラム（クリックで開く）

import math 
import numpy as np
import random
import matplotlib.pyplot as plt
import scipy.stats as st
from scipy.stats import norm


dt = 0.01
T = 1
sigma = 10
N = 10**5

t = 0
x = 0
sq_dt = math.sqrt(dt)

save_x = np.empty([0]) #1次元から配列生成:保存用
save_x2 = np.empty([0]) #1次元から配列生成:保存用
save_t = np.empty([0])

for i in range(N):
    x = 0
    t = 0
    while True:
        
        dx = random.gauss(0,sigma*sq_dt)  # 平均、標準偏差
        x = x + dx    
        
        if i == 0:
            save_x2 = np.append(save_x2,x)
            save_t = np.append(save_t,t)
        
        if t > T:
            break
        else:
            t = t + dt
            
    if i % 1e+4 == 0:
        print(N)
            
    save_x = np.append(save_x,x)

# ヒストグラム
fig = plt.figure(figsize = (10,8))

ax = fig.add_subplot(2, 1, 1)
ax.plot(save_t,save_x2)

#xx = np.linspace(0, 2.5,501)
#ax.plot(xx, st.gamma(11, 0, 1/13.).pdf(xx))

ax = fig.add_subplot(2, 1, 2)
ax.hist(save_x, bins=100, histtype='barstacked', ec='black',density=True)

a=40
xx = np.linspace(-1*a, a, 501)
p = norm.pdf(x=xx, loc=0, scale=sigma) #平均・標準偏差
ax.plot(xx, p,lw = 2)



plt.show()

ρの推定

先程の、$\rho$の推定の式をシミュレーションに適用してみよう

\hat{\rho} = \frac{\sum x_{t-1} x_t}{\sum x_t^2 }

この分布数学的に求めるのは大変難しいため、今回はプログラムでシミュレートすることで$\hat{\rho}$の分布を求めてみる。先程は以下の式に従ったがどうなるのか・・・

\rho - \hat{\rho} \sim  N(0,\frac{\sum e_t^2}{\sum y_{t-1}^2 } )

まず、分母

\sum x_t^2

Var[\epsilon_1 + \epsilon_2] = Var[\epsilon_1] + Var[\epsilon_2] + Cov[\epsilon_1,\epsilon_2]

2つは無相関であるので、

Var[\epsilon_1 + \epsilon_2] = Var[\epsilon_1] + Var[\epsilon_2]= \sigma^2 2dt

となる。

x_2^2+ x_1^2 = 2 x_1^2 + \epsilon_2^2 = 2\epsilon_1^2 + \epsilon_2^2

一般化すれば

\sum x_t^2 = \sum i \epsilon_i^2

カイ二乗分布に近いけど違う分布になります。

分子


\sum x_{t-1} x_t

あまり変わりませんね。

$\rho$の分布

average :  0.973625348374121
median :  0.9805936049706372

平均も中央値もほぼ1になることがわかった。
また、今回$\rho = 1$の過程に対して$\hat{\rho}$の推定をしたのにも関わらず平均値は$1$にならなかった。
つまり、完璧な単位根過程であっても推定値は$1$にならないのである。

そしてこの、確率密度分布をDF分布と呼ぶ。
ヒストグラムの分解能を上げれば

これがDF分布である。正式名称Dickey-Fuller分布

シミュレーションのパラメータを変えた場合もいくつか示すと

$dt = 0.007$

\epsilon_t  \sim N(0,1)

\epsilon_t  \sim N(0,100)

\epsilon_t  \sim N(0,10000)

$dt = 0.001$に下げる

\epsilon_t  \sim N(0,1)

$dt$を下げることでサンプリング回数が増えるため、推定精度が向上するため分布の形は変わらないが、横軸が1に近づく。

以上より、推定誤差分散で正規化したものがDF分布になる。

推定誤差の分布

\rho - \hat{\rho} \sim  N(0,\frac{\sum e_t^2}{\sum x_{t-1}^2 } )

であった。
$\rho = 1$の過程に対してシミュレーション上で何度も推定を行い$\hat{\rho} $を求め、$e_t$を求めて

e_t = x_t - \hat{\rho} x_{t-1}

以下の誤差推定値の分布を求める

\hat{\sigma^2 } = \frac{1}{T} \sum e_t^2

今回のシミュレーションでは

x_{t} =  x_{t-1} + \epsilon_t

\epsilon_t \in N(0,100 )

としてシミュレーションした。ちゃんと真値$100$周辺に分布している。

推定誤差を推定誤差の推定分散で正規化した値

\tau = \frac{ \rho - \hat{\rho} }{ \frac{\sum e_t^2}{\sum x_{t-1}^2 } }

の分布は

となる。

全体のプログラム（クリックで開く）

# -*- coding: utf-8 -*-
"""
Created on Wed Oct  4 00:31:58 2023

単位根過程をシミュレート

シグマの中身を計算：分布導出

OLS推定量を計算：分布導出、正規分布と比較
"""
import math 
import numpy as np
import random
import matplotlib.pyplot as plt
import scipy.stats as st
from scipy.stats import norm
from scipy.stats import chi2

dt = 0.007
T = 1
sigma = 10
N = 10**4

t = 0
x = 0
sq_dt = math.sqrt(dt)

save_x = np.empty([0]) #1次元から配列生成:保存用
save_t = np.empty([0])

save_x1 = np.empty([0])

x2 = 0
save_x2 = np.empty([0])

x3 = 0
save_x3 = np.empty([0])

x4 = 0
save_x4 = np.empty([0])

save_x5 = np.empty([0])

save_x6 = np.empty([0])

for i in range(N):
    x = 0
    t = 0
    x2 = 0
    x3 = 0
    x4 = 0
    save_x = np.empty([0]) #1次元から配列生成:保存用
    save_t = np.empty([0])
    
    while True:
        
        dx = random.gauss(0.0,sigma*sq_dt)  # 平均、標準偏差
        
        x3 = x3 + x*(x + dx)
        
        x = x + dx 
        
        x2 = x2 + x**2 
        
        #if i == 0:
        save_x = np.append(save_x,x)
        save_t = np.append(save_t,t)
        
        if t > T:
            break
        else:
            t = t + dt
            
    if i % 1e+3 == 0:
        print(i)
            
    save_x1 = np.append(save_x1,x)
    save_x2 = np.append(save_x2,x2)
    save_x3 = np.append(save_x3,x3)
    
    x4 = x3/x2
    save_x4 = np.append(save_x4,x4)
    
    e = np.append(0,save_x)[:-1] - x4 * save_x
    x5 = sum(e**2) / T
    
    save_x5 = np.append(save_x5,x5)
    
    
    x6 = (1-x4)/math.sqrt(x5)
    save_x6 = np.append(save_x6,x6)

    
    

# 代表グラフ

fig = plt.figure(figsize=(9, 7))

ax = fig.add_subplot(1, 1, 1)
ax.hist(save_x6, bins=100, histtype='barstacked', ec='black',density=True)

ax.set_xlabel("tau")
ax.set_ylabel("f")

print("average : ", np.mean(save_x6))
print("median : ", np.median(save_x6))
#ax.set_xlim([0,np.max(save_x3)])

plt.show() 


# ヒストグラム
fig = plt.figure(figsize=(7, 20))

n=7

ax = fig.add_subplot(n, 1, 1)
ax.plot(save_t,save_x)


ax = fig.add_subplot(n, 1, 2)
ax.hist(save_x1, bins=100, histtype='barstacked', ec='black',density=True)

a=sigma*3*math.sqrt(T)
xx = np.linspace(-1*a, a, 501)
p = norm.pdf(x=xx, loc=0, scale=sigma*math.sqrt(T)) #平均・標準偏差
ax.plot(xx, p,lw = 2)


ax = fig.add_subplot(n, 1, 3)
ax.hist(save_x2, bins=300, histtype='barstacked', ec='black',density=True)

ax.set_xlabel("x^2")
ax.set_ylabel("f")
ax.set_xlim([0,np.max(save_x2)*0.5])

# a=np.max(save_x2)
# xx = np.linspace(-1*a*0.01, a, 501)
# p = chi2.pdf(x=xx/(sigma**2 * sq_dt), df = N) / sigma**2 / dt
# ax.plot(xx, p,lw = 2)


ax = fig.add_subplot(n, 1, 4)
ax.hist(save_x3, bins=100, histtype='barstacked', ec='black',density=True)
ax.set_xlabel("x(t-1) * x(t)")
ax.set_ylabel("f")

ax = fig.add_subplot(n, 1, 5)
ax.hist(save_x4, bins=600, histtype='barstacked', ec='black',density=True)

ax.set_xlabel("ρ")
ax.set_ylabel("f")

print("average : ", np.mean(save_x4))
print("median : ", np.median(save_x4))

ax = fig.add_subplot(n, 1, 6)
ax.hist(save_x5, bins=100, histtype='barstacked', ec='black',density=True)
ax.set_xlabel("sum(e^2) * 1/T")
ax.set_ylabel("f")

ax = fig.add_subplot(n, 1, 7)
ax.hist(save_x6, bins=100, histtype='barstacked', ec='black',density=True)
ax.set_xlabel("tau")
ax.set_ylabel("f")


plt.show()

トレンドのあるランダムウォーク

x_{t} =  x_{t-1} + 0.05 + \epsilon_t

\epsilon_t \in N(0,100 )

としてシミュレーション場合

線形トレンドが見られる。

ランダムウォークではなくトレンドのあるARモデルでも同じような時系列モデルが作れる

x_t = 0.9 x_{t-1}  + 0.05t + \epsilon_t

x_t = 0.98 x_{t-1}  + 0.05t + \epsilon_t

x_t = 0.7 x_{t-1}  + 0.05t + \epsilon_t

ARモデル＋線形トレンドの場合、分散は時刻によらなくなるので、各線のばらつき具合は時刻$t$によらなくなっている。

$\rho$の計算結果

x_t =  x_{t-1}  + 0.05 + \epsilon_t

\epsilon_t \in N(0,100 ),T=10^4

このトレンドのあるブラウン運動に対して

\hat{\rho} = \frac{\sum x_{t-1} x_t}{\sum x_t^2 }

トレンドのないブラウン運動に対する係数の推定値を計算して、分布を求める

\hat{\rho} = \frac{\sum x_{t-1} x_t}{\sum x_t^2 }

x_t =  x_{t-1}  + 0.3 + \epsilon_t

\epsilon_t \in N(0,100 ),T=10^4

これまでの方法で$\rho$を計算するのは問題があるといえる。

OLS推定量を使用する方法。

x_t = \rho x_{t-1} + \beta t  + \epsilon_t

この過程に対して

S = \sum ( x_t - \hat{x_t})^2

を最小化するように$\rho,\beta$を推定する。
推定値の分布を求める

x_t =  0.3  x_{t-1} + 0.5 t  + \epsilon_t

$0.3$,$0.5$がちゃんと推定できている

x_t =  x_{t-1} + 0.05  + \epsilon_t

先程のDF分布とは全く異なる分布となった。

自力で考えてもわからんなぁ～

ADF,PPについて

単位根過程の中でも一番簡単な過程のみにしかDF検定は使えない。
DF検定を拡張したものがADF検定である。

今後書くかも

RでDF分布を出してみる

今後書くかも

まとめ

DF分布をシミュレーションで実際に求めることができた。

これがDF分布である。正式名称Dickey-Fuller分布