時系列データの解析を二つのデータに対して実施します。
基本的には、RawDataの確認、トレンドと季節性の排除、モデルのフィッティングの順にやっていきます。
データ①simulated data
Rで生成したデータについて分析してみます。
データを作成するコードは下記。
データのプロットと、自己相関関数を併せて出しておきます。
#### Generate time seriese data with a linear trend
>time <- 1:200
>data.ar <- arima.sim(model=list(ar=c(0.6, 0.3)), n=200, sd=1)
>data <- data.ar + 30 + 0.05 * time
>par(mfrow=c(2,1))
>plot(time, data, type="l", main="Raw data plot")
>acf(data, main="Sample autocorrelation function")
```r
出てくるデータはこんな感じです。
グラフとコレログラム(Correlogram)から線形トレンドがありそうなので、こいつを取り除きます。
```r
## Remove the trend
>linear.model <- lm(data~time)
>summary(linear.model)
Call:
lm(formula = data ~ time)
Residuals:
Min 1Q Median 3Q Max
-3.7143 -1.0238 -0.0246 1.0379 3.7570
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 27.098000 0.207381 130.7 <2e-16 ***
time 0.065674 0.001789 36.7 <2e-16 ***
>residual.series <- data - linear.model$fitted.values
>par(mfrow=c(3,1))
>plot(residual.series, main="Residual series")
>acf(residual.series, main="ACF of Residual series")
>pacf(residual.series, main="PACF of Residual series")
時系列プロットから定常性が確認でき偏自己相関から、AR(2)が妥当なオーダーだと判断出来ますので、早速モデル化。
## Fit an AR(2) model to the data
model.ar <- arima(residual.series, order=c(2,0,0))
par(mfrow=c(3,1))
plot(model.ar$residuals, main="Residual series")
acf(model.ar$residuals, main="ACF of Residual series")
pacf(model.ar$residuals, main="PACF of Residual series")
ARIMA()関数は残差をプロットしてくれます。
プロットを見る限り、残差に相関はなさそうなのでモデルが妥当だったということでしょう。
データ②Temperature Data
湖の温度を測ったデータを使っての時系列分析です。
データは私のGitHubに置いてあります。
GitHub:
[https://github.com/tkazusa/Sample_code/tree/master/Data]
> #Data import
> tempeature <- read.csv("temperature.csv", header=TRUE)
> t <- 1:nrow(tempeature)
> x <- tempeature$Temperature
> par(mfrow=c(2,1))
> plot(t,x, type="l", main="Raw data plot")
> acf(x,main="Temperature data ACF")
図から、先ほどの例とは異なり季節性が出ていることがわかります。
この季節性を取り除くために、下記のモデルを使います。
[tex:
s_t = \beta_0 + \beta_1 sin(\frac{2 \pi t}{12}) + \beta_2 cos(\frac{2 \pi t}{12})
]
> model.season <- lm(x~sin(2*pi*t/365) + cos(2*pi*t/365))
> summary(model.season)
Call:
lm(formula = x ~ sin(2 * pi * t/365) + cos(2 * pi * t/365))
Residuals:
Min 1Q Median 3Q Max
-10.3755 -2.1802 -0.0533 2.4349 9.0749
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.4263 0.1029 62.47 <2e-16 ***
sin(2 * pi * t/365) -2.5477 0.1455 -17.51 <2e-16 ***
cos(2 * pi * t/365) -3.7933 0.1455 -26.07 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.404 on 1092 degrees of freedom
Multiple R-squared: 0.4746, Adjusted R-squared: 0.4736
F-statistic: 493.2 on 2 and 1092 DF, p-value: < 2.2e-16
> residual.series <- x - model.season$fitted.values
> par(mfrow=c(3,1))
> plot(residual.series, main="Residual series", type="l")
> acf(residual.series, main="ACF of Residual series")
> pacf(residual.series, main="PACF of Residual series")
プロットからResidualデータは定常性があり、PACFからAR(1)モデルが妥当だと判断できますので、モデル化します。
>## Fit an AR(1) model
> model.ar <- arima(residual.series, order=c(1,0,0))
> model.ar
Call:
arima(x = residual.series, order = c(1, 0, 0))
Coefficients:
ar1 intercept
0.4861 0.0026
s.e. 0.0264 0.1746
sigma^2 estimated as 8.829: log likelihood = -2746.36, aic = 5498.72
> par(mfrow=c(3,1))
> plot(model.ar$residuals, main="Residual series")
> acf(model.ar$residuals, main="ACF of Residual series")
> pacf(model.ar$residuals, main="PACF of Residual series")
残差プロットやPACFより綺麗にモデル化できていると判断できます。