suzuzusu日記

(´・ω・`)

幾何ブラウン運動をARIMAで予測

忘備録としてまとめておく。

幾何ブラウン運動

以下のような確率微分方程式(SDE)によって定義されるものを幾何ブラウン運動と言います。金融工学のモデルとしてよく用いられます。

d S_t = \mu S_t dt + \sigma S_t dB_t

初期値をS_0とすると解は次のようになることが知られています。

 S_t = S_0 \exp( (\mu - \frac{\sigma ^2}{2})t + \sigma B_t) )

シミュレーションすると次のようになります。

# geometric Brownian motion
T = 1.0
N = 10000
dt = T/N
sigma = 2
mu = 1
X0 = 1
t = seq(from = 0.0, to = T, length.out = N)
dw = sqrt(dt) * rnorm(N)
dw[1] = 0.0
W = cumsum(dw)
X = X0 * exp((sigma - 0.5 * mu^2) * t + (mu*W))
plot(t, X, type="l")

f:id:suzuzusu:20191124020747p:plain
幾何ブラウン運動

ARIMA

ARIMAはRのforecastパッケージの auto.arima を使います。

library(forecast)

# forecast step
h = 1000

# ARIMA
model <- auto.arima(X[1:(N-h)])
summary(model)

ARIMA(2, 1, 2)が選択されました。以下がsummaryの出力です。

ARIMA(2,1,2) 

Coefficients:
         ar1      ar2      ma1     ma2
      0.7916  -0.7906  -0.8089  0.8230
s.e.  0.0604   0.0804   0.0559  0.0746

sigma^2 estimated as 0.001461:  log likelihood=16609.56
AIC=-33209.12   AICc=-33209.11   BIC=-33173.59

Training set error measures:
                       ME       RMSE        MAE       MPE      MAPE      MASE
Training set 0.0005448627 0.03820814 0.02565005 0.0147103 0.7971145 0.9991512
                     ACF1
Training set 0.0009177808

結果

予測の結果を次に示します。

# small plot
plot(forecast(model, h=h), xlim=c(0,N), ylim=c(0,15))
par(new=T)
plot(c((N-h):N), X[(N-h):N], xlim=c(0,N), ylim=c(0,15), ann=F, type='l', col = "red")
legend("topleft", legend = c("input", "true", "forecast"), lty=c(1, 1, 1), col = c("black", "red", "blue"))

f:id:suzuzusu:20191124021654p:plain

こちらが拡大図です。

# large plot
plot(forecast(model, h=h), xlim=c(N-h-200,N), ylim=c(0,15))
par(new=T)
plot(c((N-h):N), X[(N-h):N], xlim=c(N-h-200,N), ylim=c(0,15), ann=F, type='l', col = "red")
legend("topleft", legend = c("input", "true", "forecast"), lty=c(1, 1, 1), col = c("black", "red", "blue"))

f:id:suzuzusu:20191124021711p:plain
拡大図

実装まとめ

library(forecast)

set.seed(0)

# geometric Brownian motion
T = 1.0
N = 10000
dt = T/N
sigma = 2
mu = 1
X0 = 1
t = seq(from = 0.0, to = T, length.out = N)
dw = sqrt(dt) * rnorm(N)
dw[1] = 0.0
W = cumsum(dw)
X = X0 * exp((sigma - 0.5 * mu^2) * t + (mu*W))
plot(t, X, type="l")

# forecast step
h = 1000

# ARIMA
model <- auto.arima(X[1:(N-h)])

# small plot
plot(forecast(model, h=h), xlim=c(0,N), ylim=c(0,15))
par(new=T)
plot(c((N-h):N), X[(N-h):N], xlim=c(0,N), ylim=c(0,15), ann=F, type='l', col = "red")
legend("topleft", legend = c("input", "true", "forecast"), lty=c(1, 1, 1), col = c("black", "red", "blue"))

# large plot
plot(forecast(model, h=h), xlim=c(N-h-200,N), ylim=c(0,15))
par(new=T)
plot(c((N-h):N), X[(N-h):N], xlim=c(N-h-200,N), ylim=c(0,15), ann=F, type='l', col = "red")
legend("topleft", legend = c("input", "true", "forecast"), lty=c(1, 1, 1), col = c("black", "red", "blue"))

参考