時系列解析① ARモデル

時系列データの解析を二つのデータに対して実施します。
基本的には、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")

出てくるデータはこんな感じです。
グラフとコレログラム(Correlogram)から線形トレンドがありそうなので、こいつを取り除きます。

## 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")

図から、先ほどの例とは異なり季節性が出ていることがわかります。

この季節性を取り除くために、下記のモデルを使います。

> 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より綺麗にモデル化できていると判断できます。