2018-07-27

【時系列】ARIMA（Python）

沖本本時系列 Python

時系列データをPythonで扱うのって難しく感じるのは私だけですかねえ。

ARIMAモデル

statsmodelsのarima_modelでできます。データは2015年1月から2018年7月までの日経225データです。

import pandas as pd
import numpy as np
from statsmodels.tsa import arima_model
import matplotlib.pyplot as plt

#データ
nikkei225 = pd.read_csv('C:\\Users\\USER\\Documents\\R\\nikkei225.csv', engine='python')
data = np.array(nikkei225['終値'])

#ARIMA
results = arima_model.ARIMA(data,order = [1,1,1]).fit()

#図と予測
plt.clf()
plt.plot(data)
plt.plot(results.predict(start=0,end=50))
plt.legend(['data','predicted'])

f:id:imakoto0323:20180709134616p:plain

参考

qiita.com

2018-07-25

【時系列】ARIMA（R）

時系列沖本本 R

前回のARIMAモデルをRで実装していきたいと思います。

使用するデータ

2015年1月から2018年7月までの日経225の終値を利用したいと思います。データの作成方法は下記のサイトを参考にしました。

gist.github.com

ちなみにPythonやらで使用したいためcsvでいったん出力しています。

ARIMAモデルを入れてみる

今回は上記の終値を対数をとっています。

#data
nikkei225 <- read.csv("C:/Users/USER/Documents/R/nikkei225.csv",row.names=1)

#ARIMAモデル 
#install.packages("forecast")
library(forecast)
Arima <- auto.arima( log(nikkei225$終値),
                     ic="aic",
                     stepwise=F,
                     approximation=F,
                     max.p=5, 
                     max.q=5)

Arima
#Series: log(nikkei225$終値) 
#ARIMA(1,1,1) 
#
#Coefficients:
#  ar1      ma1
#0.6933  -0.7464
#s.e.  0.2040   0.1890
#
#sigma^2 estimated as 0.0001681:  log likelihood=2518.12
#AIC=-5030.24   AICc=-5030.21   BIC=-5015.97

#予測
yosoku <- forecast(Arima,level=c(50,95),h=30)
yosoku
#Point Forecast    Lo 50     Hi 50    Lo 95    Hi 95
#862       9.989377 9.980632  9.998122 9.963966 10.01479
#863       9.989555 9.977511 10.001598 9.954559 10.02455
#864       9.989678 9.975242 10.004113 9.947731 10.03162
#865       9.989763 9.973386 10.006141 9.942173 10.03735
#866       9.989822 9.971775 10.007869 9.937381 10.04226
#・・・

#図示
plot(forecast(Arima,level=c(50,95),h=30))
# 信頼区間を50%, 95%で設定して表示（前者の方が狭い）
# 向こう30区間の予測値と信頼区間を表示

f:id:imakoto0323:20180709132333p:plain

SARIMAもやってみる

と思ったのですが、forecastパッケージのauto.arimaが勝手に季節調整してくれる？みたいなので割愛します。

ARIMAX

auto.arimaの変数にxregを入れることで説明変数を増やすことができます。

#ARIMAX
Arimax <- auto.arima( log(nikkei225$終値),
                     xreg=nikkei225$Friday,
                     ic="aic",
                     stepwise=T,
                     approximation=T,
                     )

参考

tjo.hatenablog.com

logics-of-blue.com

2018-07-24

【時系列】ARIMA（理論編）

沖本本では2章でARMAまで説明しているのにARIMAは第5章まで待たないといけません。

本当は単位根とか説明しないといけないのかもしれませんが先に記事にしたいと思います。というのもこのARIMAモデルはビジネスの場面で本当によく見かけるからです。

今回は以下の本も参考にしてます。

時系列データの基本構造

時系列データ＝短期の自己相関＋周期的変動＋トレンド＋外因性＋ホワイトノイズという基本構造で時系列データはできています。これらの問題に周期的変動・・・ARMA成分周期的変動・・・季節成分トレンド・・・差分をとることで消す外因性・・・ARIMAXモデルで対応できます。

ARMAモデル

ARMA（自己回帰移動平均）はARモデルとMAモデルを組み合わせたモデルで自己回帰を柔軟に表現できるものでした。

f:id:imakoto0323:20180709124000p:plain

ただし、ARMAモデルは非定常過程には用いることができません。

和分過程（単位根過程）

原系列が非定常過程でも差分過程は定常過程になる場合があります。そのような過程を和分過程（単位根過程）といいます。

d期先との差分をとる和分過程をd階の和分過程といいます。

ARIMA（自己回帰和分移動平均）モデル

ARMAモデルでは非定常過程に対して利用できません。しかし和分過程に対しては差分をとって定常過程に変換してARMAモデルを利用できます。

差分をとってARMAモデルを利用するモデルをARIMAモデルといいいます。 p次先のARモデル、d階の和分過程、q次のMAモデルを合わせたARIMAモデルをARIMA(p,d,q)と表します。

SARIMA*1モデル

通常のARIMAモデルはとなりまたは数期前のデータと差分をとります。しかし季節変動がある場合にその影響を考えるために季節階差をとります。これをSARIMA（季節性ARIMA）モデルといいます。

例えば12か月周期で考えた場合、「先月も高ければ今月も高いだろう」と通常のARIMAモデルでは前月と当月との自己相関をモデル化します。 f:id:imakoto0323:20180709163701p:plain

しかしSARIMAモデルでは「1月は毎年高いだろう」と毎年同じ月の自己相関をモデル化します。 f:id:imakoto0323:20180709164136p:plain

ARIMAXモデル

ARIMAモデルに外生変数を加えた形のモデルをARIMAXモデルといいます。外生変数にはイベントや曜日等を入れることができます。

*1:「えすありま」と読む人も「さりま」と読む人もしっているのですがどちらが正解なんでしょうか

2018-07-22

【時系列】MA・AR・ARMA（理論編とR）

時系列有斐閣「統計学」統計学 R 沖本本

沖本本第２章のお話です。

有斐閣「統計学」の第12章にも載っています。

さらに今回はこちらも参考にいたしました。

はじめに

正直このあたりから少しずつ難しくなってくるころだと思います。最初読んだときはほとんど意味がわかりませんでした。

今回は要は自己相関をどうモデリングするかがテーマです。その方法は2種類あって1つ目は同じ変数（ここではb）を使って

f:id:imakoto0323:20180706120758p:plain

というあらわし方。こうすると、共通のｂによって $y_t$ と $y_{t-1}$ が相関を持つようにできます。

もうひとつの方法が $y_t$ を $y_{t-1}$ を含んだ式にすること。

f:id:imakoto0323:20180706121028p:plain

こうすることで $y_t$ と $y_{t-1}$ は相関を持たせることができます。

1つ目の考え方でモデリングしたのがMA過程で2つ目の考え方でモデリングしたのがAR過程です。

MA（移動平均）過程

f:id:imakoto0323:20180706120817p:plain

の考えでモデリングします。ホワイトノイズを拡張したモデリングとなります。 1次のAR過程をMA(1)過程とすると、

f:id:imakoto0323:20180706120916p:plain

ちなみに $y_{t-1}$ は

f:id:imakoto0323:20180706120952p:plain

となり、共通の $\epsilon_{t-1}$ を持つことになります。

MA過程は自己相関をもつので、 $\Theta_1$ が正の場合、1次の正の自己相関が生まれ、その $\Theta_1$ の値が1に近づくほど強くなるのでグラフが滑らかになります。一方、 $\Theta_1$ が負の場合、負の自己相関が生じその $\Theta_1$ 値が－1に近づくほどよりぎざぎざしたグラフとなります。

MA過程のRでの実装

ma_1 <- arima.sim(n=200,list(ma=1)) #MA次数は1、係数も1
plot(ma_1)
ma_2 <- arima.sim(n=200,list(ma=-1)) #MA次数は1、係数は-1
plot(ma_2)

f:id:imakoto0323:20180706122313p:plain f:id:imakoto0323:20180706122410p:plain

$\Theta_1$ ＞0（上）のほうがグラフが滑らかで $\Theta_1$ ＜0（下）のほうはギザギザしています。

AR（自己回帰）過程

f:id:imakoto0323:20180706121039p:plain

をホワイトノイズをもちいてモデリングすると、

となります。ちなみに過去の情報をもとに確定的に定まるのは f:id:imakoto0323:20180706121441p:plain 過去の情報とは無関係系に確率的に新たな情報を生む部分 f:id:imakoto0323:20180706121610p:plain ということで $\epsilon_1$ のみが新しい情報となります。

AR過程のRでの実装

ar_1<-arima.sim(n=200,list(ar=0.5)) # AR次数は1、係数は0.5
plot(ar_1)
ar_2<-arima.sim(n=200,list(ar=-0.5)) # AR次数は1、係数は-0.5
plot(ar_2)
ar_3<-arima.sim(n=200,list(ar=1)) # AR次数は1、係数を1にするとエラーとなる
#Error in arima.sim(n = 200, list(ar = 1)) : 
#  モデルの 'ar' 部分が定常ではありません 
'''
[f:id:imakoto0323:20180706123013p:plain]
[f:id:imakoto0323:20180706123040p:plain]
係数[tex:\mid \phi_1 v]が1を超えると、過程が指数的に上昇し**爆発的**と呼ばれる状態になります。[tex:\mid \phi_1 \mid]が1未満の場合定常となります。

ARMA（自己回帰移動平均）過程

要はAR過程＋MA過程です。 ARMA（p,q）過程は

f:id:imakoto0323:20180706123956p:plain

経済分析では系列の変動を説明するのはARの部分であり、MAは付属部分となります。

ARMAモデルの推定はOLS（最小二乗法）か最尤法で推定します。

2018-07-18

【時系列】自己相関（Python）

沖本本時系列統計学有斐閣「統計学」 Python

自己相関だけです。あとはわからなかった・・・。

自己相関

import pandas as pd
import numpy as np
from numpy.random import *
import matplotlib.pyplot as plt
import statsmodels.api as sm

#データの読込
data = pd.read_csv("C:\\Users\\imoto-mk\\Documents\\meimoku_GDP.csv")
#グラフを描く
plt.plot(data["GDP"])

#自己相関係数
acf = sm.tsa.stattools.acf(data["GDP"], nlags=40)
print(acf)
#[ 1.00000000e+00  6.79758886e-01  2.94915245e-01 -3.45919123e-02
# -2.85414892e-01 -4.40794724e-01 -4.20337187e-01 -2.47093838e-01
# -1.24956591e-01  4.60227321e-02  1.28837778e-01 -2.38933734e-02
# -1.74265306e-01 -2.14610097e-01 -1.28310049e-01 -3.37206520e-02
#  4.23469306e-02  1.56290730e-01  1.71033660e-01  1.65181535e-01
#  1.17537320e-01 -6.54451119e-05 -7.93371218e-02 -9.45336256e-02]

f:id:imakoto0323:20180626152116p:plain

2018-07-17

【時系列】自己共分散と自己相関（R）

沖本本時系列統計学有斐閣「統計学」 R

前回の実装です。自己共分散と自己相関の実装となります。

Rで実装

データは1994年から2017年の名目GDPです。

#データ読込
df <- read.csv("C:\\Users\\imoto-mk\\Documents\\meimoku_GDP.csv")
df
#   Year      GDP
#1  1994 502636.2
#2  1995 516406.5
#3  1996 528766.4
#4  1997 533338.2
#5  1998 526013.4
#6  1999 521988.3
#7  2000 528512.7
#8  2001 519073.5
#9  2002 514764.4
#10 2003 517930.6
#11 2004 521180.2
#12 2005 525692.2
#13 2006 529076.6
#14 2007 530997.3
#15 2008 509465.8
#16 2009 492070.4
#17 2010 499281.0
#18 2011 494017.2
#19 2012 494478.0
#20 2013 507246.0
#21 2014 518468.5
#22 2015 533897.3
#23 2016 539351.1
#24 2017 548696.1

#プロットしてみる
ts.plot(df, ylim=c(400000, 600000))

f:id:imakoto0323:20180626141328p:plain

自己相関

acf関数で実装でき、引数 type には "correlation" (自己相関), "covariance" (自己共分散), "partial" (偏相関) から選択できるようになっています。

#自己相関
acf(df[,2], type = "correlation")

f:id:imakoto0323:20180626142045p:plain

自己相関はなさそうですね。

自己共分散

#自己共分散
acf(df[,2], type = "covariance")

f:id:imakoto0323:20180626142305p:plain

自己相関検定

#自己相関検定(Ljung-Box検定)
#帰無仮説は「自己相関関係ではない」です。
Box.test(df[,2],lag=1,type="L")

f:id:imakoto0323:20180626143027p:plain

自己相関はなさそうです。

2018-07-16

【時系列】時系列の基本統計量（理論編）

有斐閣「統計学」統計学時系列沖本本

沖本本P6の内容です。

森棟本有斐閣「統計学」第12章でも扱われております。

自己共分散

自己共分散は、同一の時系列データの別の２時点のデータ間の共分散のことです。１期間離れたデータを１次のデータとすると、１次の共分散は f:id:imakoto0323:20180626155425p:plain となります。

自己共分散の性質

自己共分散は共分散と同じく１次の共分散が正の場合（期待値を基準に）同じ方向に動き、負の場合は逆に動きます。また１次の共分散がばらばらに動く（規則性はない）といえます。

ｋ次の自己共分散

f:id:imakoto0323:20180626160434p:plain

これをｋの関数としてみると自己共分散関数と言われ、自己共分散関数は正定値になります。正定値の解説はこちら。

datahotel.io

自己相関係数

自己共分散を、値によって変化しないようにしたものが自己相関係数です。

f:id:imakoto0323:20180626161407p:plain

自己相関係数は単に自己相関とよくいいます。自己相関をグラフに書いたものをコレログラムといい、下のような図です。

f:id:imakoto0323:20180626145438p:plain

確率過程

時系列データでは本来連続のはずで、観測点は無限にあります。しかしながら観測できるのは１点でそこから予測をするのは困難を極めます。そこで、今観測されている時系列データをある確率変数列からの１つの実現値とみなし、この確率変数の列の生成過程に仮定を置きます。これを確率過程といいます。

自己相関の検定

データが自己相関を持っている場合、その自己相関を踏まえた時系列モデルを記述できそのモデルを予測などに用いることができます。一方で自己相関を持っていないと時系列分析でできることが限られてしまいます。

つまりデータが自己相関を持っていることが大切なのです。

よく利用されるのはかばん検定というもので帰無仮説はすべての時点の自己相関が０、対立仮説は１つでも自己相関が０でない関係の２時点がある、です。

偏自己相関

$y_t$ と $y_{t-k}$ の自己相関のうち、その間のt-1期からt-k-1期の影響を取り除いたものです。

いっかくのデータサイエンティストをいく

１からプログラミングとデータサイエンスを独習したい

【時系列】ARIMA（Python）

ARIMAモデル

参考

【時系列】ARIMA（R）

使用するデータ

ARIMAモデルを入れてみる

SARIMAもやってみる

ARIMAX

参考

【時系列】ARIMA（理論編）

時系列データの基本構造

ARMAモデル

和分過程（単位根過程）

ARIMA（自己回帰和分移動平均）モデル

SARIMA*1モデル

ARIMAXモデル

【時系列】MA・AR・ARMA（理論編とR）

はじめに

MA（移動平均）過程

MA過程のRでの実装

AR（自己回帰）過程

AR過程のRでの実装

ARMA（自己回帰移動平均）過程

【時系列】自己相関（Python）

自己相関

【時系列】自己共分散と自己相関（R）

Rで実装

自己相関

自己共分散

自己相関検定

【時系列】時系列の基本統計量（理論編）

自己共分散

自己共分散の性質

ｋ次の自己共分散

自己相関係数

確率過程

自己相関の検定

偏自己相関