幾何ブラウン運動をARIMAで予測
忘備録としてまとめておく。
幾何ブラウン運動
以下のような確率微分方程式(SDE)によって定義されるものを幾何ブラウン運動と言います。金融工学のモデルとしてよく用いられます。
初期値をとすると解は次のようになることが知られています。
シミュレーションすると次のようになります。
# 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")
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"))
こちらが拡大図です。
# 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"))
実装まとめ
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"))