1.0 Rで時系列データを扱う際の基礎

tsクラス：データに時間に関する情報を加えて構造化したもの
ts.union()：複数の時系列データを統合する関数
tsp()：開始時点、終了時点、周期を取得する関数
as.Date()：文字列を年月日を表すDateオブジェクトに変換する関数

使えそうなやつ↓

days = seq(from = as.Date('2000-01-01'), to = as.Date('2000-01-31'), by = 'day')
days

###出力###
 [1] "2000-01-01" "2000-01-02" "2000-01-03" "2000-01-04" "2000-01-05" "2000-01-06" "2000-01-07" "2000-01-08"
 [9] "2000-01-09" "2000-01-10" "2000-01-11" "2000-01-12" "2000-01-13" "2000-01-14" "2000-01-15" "2000-01-16"
[17] "2000-01-17" "2000-01-18" "2000-01-19" "2000-01-20" "2000-01-21" "2000-01-22" "2000-01-23" "2000-01-24"
[25] "2000-01-25" "2000-01-26" "2000-01-27" "2000-01-28" "2000-01-29" "2000-01-30" "2000-01-31"

1.1 ナイル川における年間流量

1.1.0 下調べ

横軸時間とヒストグラムをプロットする。
自己相関係数をプロットする。
周波数領域変換を行う。

#---
#　グラフィックパラメータの変更　
#---

oldpar = par(no.readonly = T) 
# par()でグラフィックパラメータを変更する際は、元々のパラメータを
# 適当な変数に保存しておいて、作図終了後にもとに戻しておくのが良い。

par(mfrow  = c(2,2))               # 2行2列
par(oma    = c(0, 0, 0, 0))        # 外側余白
par(mar    = c(4, 4, 2, 1))        # プロットの四隅に置かれる余白のサイズ(?)
par(family = "HiraginoSans-W4")    # 日本語

#---
#　プロット
#---

plot(Nile, main = 'plot')

#---
#　ヒストグラム
#---

hist(Nile, main = 'histgram', xlab = '流水量')

#---
#　自己相関係数
#---

acf(Nile, main = 'コレログラム')

#---
#　周波数スペクトルを描画
#---

plot.spectrum = function(dat, label = '', main = '', 
                         y_max = 1, tick = c(8, 4), unit = 1){
  #
  #周波数スペクトルを描画する関数。
  #
  #Parameters
  #
  #dat：データ
  #lab,main：何ヶ月周期かを示すラベルとグラフのタイトル
  #y_max：スペクトルの最大値(一回1でやってから小さくしていく)
  #tick：横軸で目盛を打つ時点
  #
  
  # データの周波数領域変換
  dat_FFT = abs(fft(as.vector(dat)))  # 高速フーリエ変換を行う関数fft()
  
  # グラフ横軸(周波数)の表示に関する準備
  dat_len  = length(dat_FFT)
  freq_tick = c(dat_len, tick, 2)
  
  # 周波数領域でのデータをプロット
  plot(dat_FFT / max(dat_FFT), type = 'l', main = main,
       ylab = '|規格化後の周波数スペクトル|', ylim = c(0, y_max),
       xlab = sprintf('周波数[1/%s]', label), xlim = c(1, dat_len / 2), xaxt = 'n')
  axis(side = 1, at = dat_len / freq_tick * unit + 1,
       labels = sprintf('1/%d', freq_tick), cex.axis = 0.7)
}

plot.spectrum(Nile, label = '年', main = '',
              y_max = 0.1)

par(oldpar)

データから、不規則な変動が続いていることがわかる。歴史的知見より、1899年に流水量が急減したことがわかっているが、変動に埋もれていて明確には分からない。
極端な外れ値はない。
全体的に自己相関の値は極端に小さい値をとってはいない。
特定の周期性は認められない。

1.1.1 モデルの定義、フィルタリング、評価

下調べから、Nileデータは同じような値が続くものの不規則に変動していて周期性はないことがわかっている。そこで、

データ＝レベル成分＋残差

というモデルを考え、ホルト・ウィンタース法によってフィルタリングを行う。

#---
#　フィルタリング
#---

HW_Nile = HoltWinters(Nile, 
                      beta = FALSE, # トレンド成分を考えない
                      gamma = FALSE # 周期成分を考えない
                      )

#---
#　フィルタリング値の描画
#---

mygray = '#80808080'
plot(HW_Nile, main = '', 
     col = mygray, col.predicted = 'black', lty.predicted = 'dashed')

#---
#　成分に分けて描画
#---

HW_decomp = ts.union(y = HW_Nile$x,                     # 元のデータ
                     Level = HW_Nile$fitted[,'level'],  # Level成分
                     Residuals = residuals(HW_Nile)      # 残差
                     )
plot(HW_decomp, main = '')

推定値の揺らぎがやや大きいものの概ね似たような値が続き、1899年頃から全体的な値が減少している傾向も見て取れる。また全体的に残差が大きすぎたり相関が残っている印象もない。実際に残差の自己相関をプロットしてみると、以下のようになる。

#---
#　残差の自己相関
#---

acf(residuals(HW_Nile))

1.2 大気中の二酸化炭素濃度

1.2.0 下調べ

こちらからデータを取得する。そのまま読み込むと「不正なマルチバイト文字があります」と怒られる。詳細はこちら参照。
データをts型に変換し、2014年12月までで打ち切る。
欠損値の補完
横軸時間とヒストグラムをプロットする。
自己相関係数をプロットする。
周波数領域変換を行う。

#---
#　データの読み込み
#---

Ryori = read.csv('co2_monthave_ryo.csv', fileEncoding = 'SJIS',
                 stringsAsFactors = FALSE)
y_all = ts(data = as.numeric(Ryori$二酸化炭素濃度の月平均値.綾里..ppm.),
           start = c(1987, 1), frequency = 12)
y_CO2 = window(y_all, end = c(2014, 12))

#---
#　データの読み込み
#---

Ryori = read.csv('co2_monthave_ryo.csv', fileEncoding = 'SJIS',
                 stringsAsFactors = FALSE)
y_all = ts(data = as.numeric(Ryori$二酸化炭素濃度の月平均値.綾里..ppm.),
           start = c(1987, 1), frequency = 12)
y_CO2 = window(y_all, end = c(2014, 12))

#---
#　グラフィックパラメータの変更(省略)
#---

#---
#　プロット
#---

plot(y_CO2, main = 'plot')

#---
#　ヒストグラム
#---

hist(y_CO2, main = 'histgram', xlab = '二酸化炭素濃度')

#---
#　自己相関係数
#---

acf(y_CO2, main = 'コレログラム')

#---
#　周波数スペクトルを描画
#---

plot.spectrum(y_CO2, lab = '月', main = '',
              y_max = 0.02, tick = c(12, 6))

par(oldpar)

右上がりで年周期のパターンが認められる。
極端な外れ値はない。
全体的に自己相関係数は極端に小さい値をとってはいない。
12ヶ月で周期性が認められる。

1.2.1 モデルの定義、フィルタリング、予測

下調べから、右肩上がりの傾向と年単位の周期性があるデータであることがわかっている。そこで、

データ＝レベル＋傾き＋周期

というモデルを考え、ホルト・ウィンタース法によってフィルタリングとを行う。

#---
#　フィルタリング
#---

HW_CO2 = HoltWinters(y_CO2)

#---
#　フィルタリング値の描画
#---

mygray = '#80808080'
plot(HW_CO2, main = '', 
     col = mygray, col.predicted = 'black', lty.predicted = 'dashed')

#---
#　成分に分けて描画
#---

HW_decomp = ts.union(y         = HW_CO2$x,
                     Level     = HW_CO2$fitted[,'level'] + HW_CO2$fitted[,'level'],
                     Season    = HW_CO2$fitted[,'season'],
                     Residulas = residuals(HW_CO2))
plot(HW_decomp, main = '')

つぎに、2015年の値を予測する。

#---
#　予測
#---

par(oldpar)
HW_predict = predict(HW_CO2, n.ahead = 12)
plot(HW_CO2, HW_predict, main = 'ホルト・ウィンタース法によるフィルタリングと予測',
     col = mygray, col.predicted = 'black', lty.predicted = 'dashed')

次回：基礎から分かる時系列分析 2.カルマンフィルタ(数理的な導出)