使用環境

env.jl

julia> versioninfo()
Julia Version 1.7.0
Commit 3bf9d17731 (2021-11-30 12:12 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i5-10210U CPU @ 1.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake)

(@v1.7) pkg> st Optim
      Status `C:\Users\user_name\.julia\environments\v1.7\Project.toml`
  [429524aa] Optim v1.5.0

パラメータ推定

前回はボラティリティ変動モデルとはといった話をしました。今回は、前回紹介したARCHモデルのパラメータ推定を実装していきます。

沖本本や渡部本にあるようにARCHモデルは最尤推定を使ってパラメータを推定します。ただ前回を見てわかるように、ARCHモデルはARモデルに書き換えることができます。なので、ARモデルと同じようなパラメータ推定を行って、ARCHモデルのパラメータを推定することが可能です。

今回は前回とは違い、ARモデルに書き換えることなく、直接ARCHモデルのパラメータを推定したいと思います。

方針

ARCHモデルの推定法を調べると以下の２つの方針があることがわかりました。

同時推定
二段階推定

ARCHモデルは期待値部分のモデルと分散部分のモデルを両方兼ね備えたモデルです。二段階推定とは、期待値部分のモデルパラメータを推定したあとに、分散部分のモデルパラメータを推定します。
私が調べられる範囲では、二段階推定が主流で、同時推定を行っている情報を見つけることはできませんでした。

以上のことを考慮し、実装は二段階推定で行うものとします。

同時推定について

今回の実装に当たって、PythonでARCHモデルを記述できるarchというパッケージを参考にしました。Pythonのほうは私がソースを見た限りでは二段階推定を行っています。

反対にRのrugarchは同時推定を行っているようです（「ようです」なのはソースを私が読めなかったので確認できていないという、、、）。

いろいろ調べたのですが、同時推定について言及している資料を見つけることがで来ませんでした。ただ、理屈でいえば、パラメータの更新を同時に行えばいいと思うので、

期待値部分のパラメータを初期化する
期待値部分のパラメータを用いて, 残差を計算する
分散部分のパラメータを更新
更新した分散部分のパラメータを使って期待値部分のパラメータを更新
以下2 ~ 4 を指定回数またはパラメータの変化が小さくなるまで続ける

というのが考えられるかなと。
どなたか同時推定について詳しかたがいらっしゃれば教えてください。

Pythonパッケージ arch について

archは、期待値部分に含められるモデルはあらかじめ決まっています。ざっくりいうとARモデルとそれのちょっとした拡張までです。それ以外のモデルを組み込みたいときは、statsmodelsなどのほかのパッケージを使って、まず期待値のモデル化を行い、残差を求めます。これを所与として分散部分のモデル化を学習させるというフローです。

https://pypi.org/project/arch/ より抜粋

例1）AR-ARCH
ARは、archだけで学習が完了します

例2）ARMA-ARCH
ARMAはarchでは書けないので、archだけではなく、statsmodelsも使います。

Pythonのコードが見たい方はこちらの記事をご覧ください。

対数尤度関数

以下、期待値部分のモデルを $\mu_t$, 分散部分のモデルを $u_t$ とします。また今後特に断りなく$\theta$ を使うときはパラメータをまとめたものと考えてください。

まず私たちの出発地点は、$\epsilon_t$ をホワイトノイズとした観測値 $y$ についての条件付き期待値のモデル化から始まります。

$$
y_t = \mu_t + \epsilon_t \quad \epsilon ～ N( 0, \sigma )　\quad \tag{1}
$$

前回記事より、残差の二乗に自己相関が現れることを確認しました。
そこで、$\epsilon$ をそれぞれ独立な $u$ と $z ～ N(0,1)$ に分解できると仮定します。すなわち

$$
\epsilon_t = \sqrt{u_t} z_t \quad z_t ～ N(0,1) \quad \tag{2}
$$

ここで、$u_t$ は次のような式で表されます。

$$
u_t = w + \Sigma_{p=1}^{P} \alpha_{t} \epsilon_{t-p} \quad \tag{3} \
w > 0, \alpha \geq 0
$$

(1),(2)より、

\begin{align}
y_t &= \mu_t + \epsilon_t \\
\epsilon_t &= y_t - \mu_t \\
\sqrt{u_t} z_t &= y_t - \mu_t
\end{align}

$z_t ～ N(0, 1)$ であることを考えると、

$$
\epsilon_t ～ N( \mu_t, u_t )
$$

すなわち

$$
f(\epsilon_t|\theta) = -\frac{1}{\sqrt{2\pi u_t}} \exp{(-\frac{1}{2u_t}(y_t - \mu_t^2))}
$$

となります。尤度関数は確率密度の積で表されますが、通常対数をとったものを最適化します。

$$
L(X|\theta) = -\frac{T-p-1}{2} \log{2\pi}

\frac{1}{2} \Sigma_{t=1}^{T-p-1} \log{u_{t+p}}
\frac{1}{2} \Sigma_{t=1}^{T-p-1} \frac{1}{u_{t+p}}(y_{t+p} - \mu_{t+p})^2 \quad \tag{4}
$$

(4)はプログラムでは次のように書けます。

loglikelihood.jl

sum( logpdf( Normal( μ, σ ), X )

logpdfは確率密度の値に対数を取ったものです。

観測値の分布について

本記事では何の気になしに正規分布を仮定しましたが、沖本本によれば実証実験では正規分布のほかにt分布やGED（一般化誤差分布）を用いることもあるらしいです。
詳しくは沖本本や渡部本を参考にしてください。

実装

使い方の例としてmain.jlを書きました。estimateParamsOfArchがARCHモデルのパラメータ推定を行う関数です。第一引数に、期待値部分のモデルの残差の二乗を入れ、第二引数にARCHモデルの次数を選択します。

返り値は、Optim.jlの最適化結果と対数尤度関数を返します。

※AR部分の実装は過去記事のグレンジャー因果性検定の実装計量時系列分析4.3 P81～P82にあるのですが、使いづらかったので書き直しました。詳しい実装は一番下のソースnew_ar.jlに書きました。

main.jl

# size( X ) = ( 100 )
data = diff( log.( X ) )
# ar学習
ar = Q.ar(1)
ar_res = Q.fit( ar, data )
# 残差の2乗
resid2 = Q.resid( ar_res, data ).^2

# ARCH(1)
arch_res, f = estimateParamsOfArch( resid2, 1 )

# 推定されたパラメータ
arch_res.minimizer

ARCHモデルパラメータ推定

以前までの記事では、対数尤度関数もすべて自前で書いていたのですがどうしてもパラメータが収束しなかったので、今回初めてDistributions.jlの力を借りて書いてみました。

arch.jl

using Optim

"""
    estimateParamsOfArch( r², lags ) -> Optim.MultivariateOptimizationResults{ ... 略 }, function
    
　　r²　: 残差の二乗 
　　lags: archモデルの次数
"""
function estimateParamsOfArch( r², lags )
    X = r²[end:-1:1]
    param = rand( lags + 1 )
    _X = hcat( ( X[begin + i:lags + i] for i in 1:length(X) - lags )... )

    # arch
    ut( x, w ) = w' * [ 1; x ]
    
    # 対数尤度関数
    f( x, X, _X, lags ) = sum( logpdf.( Normal.( 0.0, ut( _X, x )' ), X[begin:end - lags] ) )
    
    # 準ニュートン法を使用 
    # memo 最急降下法だと学習時間がかかりすぎた
    inner_optimizer = BFGS()
    # 最適化の様子を追えるように関数値を内部に保存する
    optim_ops = Optim.Options(
        store_trace=true,
        extended_trace=true,
    )
    
    # パラメータに非負制約を加える
    lower = zeros( lags + 1 )
    upper = [ Inf for _ in 1:lags + 1 ]

    # 最適化 optim.jlはパラメータを最小化するので, 対数尤度関数にマイナス1を掛ける
    res = optimize( 
        x -> -f( x, X, _X, lags ),
        lower,
        upper,
        rand( lags + 1 ),
        Fminbox(inner_optimizer),
        optim_ops
    )
    res, f
end

ARモデルの書き直し

下記コードのloglikelihood loglikelihood2 はうまく収束させることができなかったのですが、供養として。

new_ar.jl


module Q

using Statistics ; mean
using LinearAlgebra ; dot
using Optim
#using LineSearches
using Distributions

abstract type ConditionalMean end
abstract type ConditionalVolatility end

"""
    arch{T <: Integer} <:ConditionalVolatility

    lags of arch model
"""
struct arch{T <: Integer} <:ConditionalVolatility
    q::T
end

"""
    ar{T <: Integer} <:ConditionalMean

    lags of ar model
"""
struct ar{T <: Integer} <:ConditionalMean
    p::T
end

"""
    q         : lags
    w         : const
    α        : coef
    mean      : params of conditional mean model
    optim_res : return optimize
    ll        : likelihood value
"""
struct ArchResults
    q
    w
    α
    mean # 条件付き期待値のモデルパラメータ
    optim_res
    ll   # 尤度の値を格納する配列
end

"""
    p  : lags
    ϕ : coef
    σ : var
"""
struct ArResults
    p
    ϕ
    σ
end

"""
    μ : const
"""
struct ConstantResults
    μ
end

"""
    lenght of params
"""
Base.length( v::arch ) = v.q + 1

@doc raw"""
    likelihood(a::arch, m::ArResults, x, X, resid, init_var; ρ=1_000_000)

    a          : arch obj
    m          : optimized ar params 
    x          : params to be optimize. x = [ w, α₁, α₂ ... ]
    X          : samples
    resid2     : Squared difference between mean model and obs
    init_var   : variance of resid
    ρ         : penalty coef

# HACK \ はQiitaに張るために入れたものなので、実際に使うときは消す 
\```math
L(\theta) = - \frac{1}{2} \Sigma_{t=1}^{T} h_t - \frac{1}{2} \Sigma_{t=1}^{T} \frac{( y_t - \mu_t )^2}{h_t}
+ \frac{\rho}{w} + \Sigma_{q=1}^{Q} \frac{\rho}{\alpha_q}  
\```
    return log likelihood which have penalty function.
""" 
function loglikelihood1(a::arch, m::ArResults, x, X, resid2, init_var; ρ=1_000_000_000)
    # 残差を二乗し, 初期値を先頭に加える 
    u = vcat( [ init_var for _ in 1:a.q ], resid2 )
    # 条件付き期待値のモデルの初期値を設定する 
    # 沖本本では適当な定数を置くとしている
    # 初期値を0.0と設定
    X = vcat( [ 0.0 for _ in 1:m.p ], X[begin:end-m.p] )
    # 対数尤度の初期化
    ll = 0.0

    @show x
    # 尤度の計算
    for i in 1:length(u) - a.q + 1
        # 条件付き分散
        h = x[1] + x[2:end]' * u[i:a.q+i-1][end:-1:begin]
        # 定数項と第2項
        ll += ( log(2pi) + log(h) ) / 2
        # 第3項
        ll += ( X[i] - m.ϕ[1] - m.ϕ[2:end]' * X[i:m.p+i-1][end:-1:begin] )^2 / 2h 
    end
    # ペナルティ関数
    pen = ρ * ( 1 / x[1] + sum( 1 / x[i] for i=1:a.q ) )
    
    # 尤度関数値
   @show - ll#  + pen
end

function loglikelihood2(a::arch, m::ArResults, x, X, resid2, init_var)
    linear( coef, data ) = coef' * vcat( ones(), data )
    # 初期値を加える
    u = vcat( [ init_var for _ in 1:a.q ], resid2 )

    ll = 0.0 
    for i in 1:length(u) - a.q + 1
        h = linear( x, u[i:i+a.q-1][end:-1:begin] )
        x_ = linear( m.ϕ, X[i:i+m.p-1][end:-1:begin] )
        ll += loglikelihood( Normal( x_, sqrt(h) ), [X[i+a.q-1]] )
    end
    ll
end

#y2   y1 u1 u0   y3  y2 y1
#y10  y9 u9 u8  y10  y9 y8

resid(m::ArResults, X) = X[begin+m.p:end] .- predict(m, X)
resid(c::ConstantResults, X) = X .- c.μ

#===== prediction =====#
"""
    predict(m::ArResults, X)

    predict in samples
"""
function predict(m::ArResults, X)
    res = zeros(length(X) - m.p)
    for i in 0:length(res)-1
        rev_data = X[begin+i:m.p+i][end:-1:1]
        
        basis_func =[ 1 ; rev_data ]
        res[i+1] = dot( m.ϕ, basis_func ) 
    end
    res
end

#===== fit =====#
"""
    fit( m::T, X ) where T <: Real

    estimate ar model params
"""
fit( m::T, X ) where T <: Real = ConstantResults(m.μ)

"""
    fit(m::CM, v::CV, X) where {CM <:ConditionalMean, CV <:ConditionalVolatility}

    estimate arch model params
    mean       : constant or ar
    volatility : arch
"""
function fit(m::CM, v::CV, X) where {CM <:ConditionalMean, CV <:ConditionalVolatility}
    # 条件付き期待値と観測値の残差の2乗を求める
    mean_params = fit( m, X )  
    u = resid( mean_params, X ).^2
    # 残差の二乗の初期値
    u_ini = var(u)
    
    inner_optimizer = GradientDescent()
    optim_ops = Optim.Options(
        store_trace=true,
        extended_trace=true,
    )
    lower = [0.0, 0.0]
    upper = [Inf, Inf]
    
    # 最適化 optim.jlはパラメータを最小化するので, 対数尤度関数にマイナス1を掛ける
    res = optimize( 
        x -> -loglikelihood2( v, mean_params, x, X, u, u_ini ),
        lower,
        upper,
        rand( length(v) ),
        Fminbox(inner_optimizer),
        optim_ops
        )

    res, mean_params,x->loglikelihood2( v, mean_params, x, X, u, u_ini )
    #ArchResults( q, w, α, mean_params, res, ll )
end

"""
    fit( m::ar, X )

    estimate ar model params
"""
function fit( m::ar, X )
    N = size(X)[1]
    
    Φ = zeros( N-m.p, m.p+1 ) 
    # 定数項
    Φ[:,1] .+= 1
    for i in 0 : m.p-1
        Φ[ :,i+2 ] += X[ m.p-i : end-1-i ]
    end

    Y = X[ m.p+1 : end ]

    coef = inv(Φ'Φ)Φ'Y
    σ² = mean( ( Y .- (coef'Φ')' ).^2 )
    
    return ArResults(m.p , coef, σ²)
end
    
end

参考文献