経営のための統計学II

時系列データの種類

時系列分析の目的は、ほとんどの場合には時系列データそのもの（原系列と呼ばれる）の性質を明らかにすることであるが、実際の解析は、原系列に何らかの変換を施した系列に対して行われることも少なくない。そのような変換を施した系列も含めて、時系列分析で扱われるデータの種類を、以下に簡単に整理する。

原系列：時系列データそのもの
対数系列：原系列の対数値
階差系列：１時点離れたデータとの階差または差分
変化率（成長率）：１時点離れたデータから変化した割合
- 離散型成長率：１時点離れたデータの比から１を引く
- 連続型成長率：対数値の差分
季節調整済み系列：原系列から季節変動を取り除いた系列
- 略して季調済み系列と呼ぶこともある。
- 季節調整自体は非常に難しい問題であり重要な問題であるが、この授業では扱わない。

時系列プロット

時系列分析では、まず時系列データをプロットして、その性質を検討する。

以下では、例として季節調整済みの鉱工業生産指数（Indices of Industrial Production：IIP）を利用する
- 経済産業省：https://www.meti.go.jp/statistics/tyo/iip/b2015_result-2.html

鉱工業生産指数データ（csvファイル）の読み込み

# IIPデータの読み込み
iip <- read.csv("iip.csv")
head(iip)

##   yyyymm  iip
## 1 197801 63.0
## 2 197802 62.4
## 3 197803 63.9
## 4 197804 63.9
## 5 197805 64.1
## 6 197806 64.4

Rにおける時系列データの扱い

Rで時系列データを扱うのに適したオブジェクトの基本は、tsクラス（tsはtime seriesの略）のオブジェクト

これは、データに時間に関する情報を加えて構造化したもの

# データをts型に変換
x <- ts(data = iip$iip, start = c(1978, 1), frequency = 12)
x

##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 1978  63.0  62.4  63.9  63.9  64.1  64.4  64.8  65.7  66.3  66.4  66.6  67.2
## 1979  67.1  67.0  67.7  67.9  69.4  69.6  70.1  71.4  69.9  71.4  72.2  72.0
## 1980  72.6  75.3  73.3  74.5  73.9  73.3  73.2  71.1  72.1  72.6  71.6  72.6
## 1981  72.7  72.4  72.7  72.7  71.6  73.5  73.9  73.7  74.9  75.4  75.4  75.2
## 1982  74.7  74.2  74.9  73.8  73.3  74.2  73.9  74.0  74.4  72.2  74.0  73.4
## 1983  73.7  73.5  74.5  74.6  74.8  75.8  75.4  77.0  78.4  77.9  78.6  79.6
## 1984  80.0  81.9  81.4  81.3  83.0  83.3  83.2  84.5  83.5  86.0  85.9  85.3
## 1985  86.1  86.1  85.3  86.6  87.6  86.1  87.4  86.5  85.7  86.6  86.5  86.4
## 1986  86.4  86.3  86.2  86.2  86.1  86.5  86.6  84.4  87.4  86.1  84.9  87.2
## 1987  86.5  86.3  87.6  86.6  85.3  88.5  89.3  89.7  91.0  92.4  92.8  93.9
## 1988  94.5  97.2  95.2  97.7  96.8  97.0  97.4  97.9  99.6  98.1 100.1 101.0
## 1989 102.1 101.3 105.0 103.2 103.1 104.6 102.3 104.3 103.7 103.0 104.0 104.4
## 1990 103.5 104.4 106.4 106.4 107.2 107.4 108.5 109.1 108.5 110.6 109.8 109.8
## 1991 110.8 111.0 109.8 110.3 111.4 108.0 110.6 109.0 108.9 108.2 108.9 107.4
## 1992 106.5 106.1 104.4 104.1 101.6 103.9 103.5 100.3 104.4 101.0  99.6  99.0
## 1993 100.8 101.9  97.8  99.6 100.9  99.1  99.6  98.7  98.7  97.4  97.3  96.5
## 1994  95.4  97.0  97.9  99.0  98.6  99.9 100.8 101.6 101.0 102.0 102.6 103.2
## 1995 100.6 102.7 103.7 105.0 102.9 103.2 101.6 102.5 101.6 102.9 103.7 105.0
## 1996 103.1 103.3 103.1 103.3 104.1 104.2 105.0 105.5 106.6 107.0 107.4 108.4
## 1997 109.8 110.1 110.8 107.9 111.0 110.2 110.9 110.3 109.4 108.6 107.2 106.2
## 1998 107.7 105.0 102.8 102.2 101.3 101.4 101.2  99.0 100.6 100.0  99.6  99.4
## 1999 100.6 100.6 102.3 100.1 102.0 101.1 101.9 102.4 103.2 103.3 104.6 104.4
## 2000 104.4 104.1 105.9 107.3 106.9 108.5 107.9 109.4 107.2 108.9 109.3 110.7
## 2001 106.0 107.2 105.4 104.3 102.2 101.0  99.4  98.2  96.2  96.1  94.5  95.5
## 2002  94.8  96.3  97.1  96.5 100.7  99.6 100.3 100.6 101.3 101.4 101.0 100.9
## 2003 101.4 101.0 101.6 100.3 101.7 100.9 101.6 100.3 103.2 105.0 104.7 104.6
## 2004 106.3 106.1 105.6 107.4 107.4 107.7 109.0 107.8 108.0 106.4 107.4 106.0
## 2005 108.4 108.2 108.6 109.1 108.4 108.7 107.8 107.9 108.9 108.4 110.1 110.3
## 2006 110.8 110.7 111.3 113.5 111.9 113.3 113.7 114.1 114.1 115.0 115.4 115.7
## 2007 114.4 115.1 115.1 114.6 115.9 116.0 116.1 119.1 117.2 119.4 117.7 118.5
## 2008 119.1 119.4 118.3 117.6 118.2 114.9 114.7 110.7 112.0 109.3 102.0  93.6
## 2009  85.3  78.0  79.0  82.5  85.5  87.1  88.3  89.6  92.6  95.0  97.0  97.8
## 2010 100.3 100.7 100.9 102.0 101.8 101.0 102.1 102.5 104.1 101.2 102.8 103.4
## 2011 103.9 104.5  87.3  89.2  95.3  99.3 100.5 102.2 101.3 103.1 100.9 102.9
## 2012 103.3 103.1 102.9 102.4 100.6  99.8  99.3  97.8  95.7  96.0  95.1  96.4
## 2013  94.8  96.5  97.7  97.7  99.3  98.2  99.8 100.0 101.0 101.2 101.8 101.8
## 2014 103.8 102.7 104.2  99.6 101.9 100.3 100.1  99.5 100.7 100.4 100.4  99.9
## 2015 102.9  99.8  99.3  99.5  99.5 100.4 100.3  98.6 100.6 100.7  99.9  98.5
## 2016 100.1  99.2  99.7  99.3  98.5  99.2  99.8 100.5 100.7 101.0 102.0 102.0
## 2017 100.9 101.6 101.5 104.1 102.3 103.3 102.5 104.0 103.0 103.3 104.2 105.8

frequencyはデータの間隔（周期）を意味しており、月次データならfrequency = 12と指定する。
- 四半期データならfrequency = 4と指定する。

原系列

# 時系列プロット：原系列
par(family = "HiraginoSans-W3") #日本語が文字化けしないようにおまじない
plot(x, xlab = '年月', ylab = '鉱工業生産指数の原系列')

1980年代は右上がりのトレンドがあるように見える。
1990年代はバブル崩壊後に上がったり下がったりしている。
2000年代は右上がりのトレンドの後に大きな落ち込み（世界金融危機）がある
2010年代は金融危機から回復（上昇）していたが、地震の影響で低下し、その後に再度上昇傾向が見られる。

対数系列

時間と共に変動幅が変わる不均一分散を緩和する目的がある。

# 時系列プロット：対数系列
par(family = "HiraginoSans-W3") #日本語が文字化けしないようにおまじない
plot(log(x), xlab = '年月', ylab = '鉱工業生産指数の対数系列')

原系列とあまり変わらない
鉱工業生産指数の原系列は単調で、対数変換は必要ないと考えられる。

階差系列

GDPの原系列にはトレンドが含まれる。
トレンドの変化を検出する場合は、階差系列を検討する
時系列$X_t$が単純なトレンド（観測時$t$の１次式と誤差項$\epsilon_t$の和） \[ X_t = \mu + t \alpha + \epsilon_t \] で表現できる場合、階差系列は次のようになる。 \[ \Delta X_t = X_t - X_{t-1} = \alpha + (\epsilon_t - \epsilon_{t-1})\]
- 平均は定数$\alpha$
階差プロットが一定の値を中心として振動するならば、原系列には１つのトレンドがあるといえる。
階差の計算には関数diffを使う。

# 階差系列
dx <- diff(x)
dx

##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 1978        -0.6   1.5   0.0   0.2   0.3   0.4   0.9   0.6   0.1   0.2   0.6
## 1979  -0.1  -0.1   0.7   0.2   1.5   0.2   0.5   1.3  -1.5   1.5   0.8  -0.2
## 1980   0.6   2.7  -2.0   1.2  -0.6  -0.6  -0.1  -2.1   1.0   0.5  -1.0   1.0
## 1981   0.1  -0.3   0.3   0.0  -1.1   1.9   0.4  -0.2   1.2   0.5   0.0  -0.2
## 1982  -0.5  -0.5   0.7  -1.1  -0.5   0.9  -0.3   0.1   0.4  -2.2   1.8  -0.6
## 1983   0.3  -0.2   1.0   0.1   0.2   1.0  -0.4   1.6   1.4  -0.5   0.7   1.0
## 1984   0.4   1.9  -0.5  -0.1   1.7   0.3  -0.1   1.3  -1.0   2.5  -0.1  -0.6
## 1985   0.8   0.0  -0.8   1.3   1.0  -1.5   1.3  -0.9  -0.8   0.9  -0.1  -0.1
## 1986   0.0  -0.1  -0.1   0.0  -0.1   0.4   0.1  -2.2   3.0  -1.3  -1.2   2.3
## 1987  -0.7  -0.2   1.3  -1.0  -1.3   3.2   0.8   0.4   1.3   1.4   0.4   1.1
## 1988   0.6   2.7  -2.0   2.5  -0.9   0.2   0.4   0.5   1.7  -1.5   2.0   0.9
## 1989   1.1  -0.8   3.7  -1.8  -0.1   1.5  -2.3   2.0  -0.6  -0.7   1.0   0.4
## 1990  -0.9   0.9   2.0   0.0   0.8   0.2   1.1   0.6  -0.6   2.1  -0.8   0.0
## 1991   1.0   0.2  -1.2   0.5   1.1  -3.4   2.6  -1.6  -0.1  -0.7   0.7  -1.5
## 1992  -0.9  -0.4  -1.7  -0.3  -2.5   2.3  -0.4  -3.2   4.1  -3.4  -1.4  -0.6
## 1993   1.8   1.1  -4.1   1.8   1.3  -1.8   0.5  -0.9   0.0  -1.3  -0.1  -0.8
## 1994  -1.1   1.6   0.9   1.1  -0.4   1.3   0.9   0.8  -0.6   1.0   0.6   0.6
## 1995  -2.6   2.1   1.0   1.3  -2.1   0.3  -1.6   0.9  -0.9   1.3   0.8   1.3
## 1996  -1.9   0.2  -0.2   0.2   0.8   0.1   0.8   0.5   1.1   0.4   0.4   1.0
## 1997   1.4   0.3   0.7  -2.9   3.1  -0.8   0.7  -0.6  -0.9  -0.8  -1.4  -1.0
## 1998   1.5  -2.7  -2.2  -0.6  -0.9   0.1  -0.2  -2.2   1.6  -0.6  -0.4  -0.2
## 1999   1.2   0.0   1.7  -2.2   1.9  -0.9   0.8   0.5   0.8   0.1   1.3  -0.2
## 2000   0.0  -0.3   1.8   1.4  -0.4   1.6  -0.6   1.5  -2.2   1.7   0.4   1.4
## 2001  -4.7   1.2  -1.8  -1.1  -2.1  -1.2  -1.6  -1.2  -2.0  -0.1  -1.6   1.0
## 2002  -0.7   1.5   0.8  -0.6   4.2  -1.1   0.7   0.3   0.7   0.1  -0.4  -0.1
## 2003   0.5  -0.4   0.6  -1.3   1.4  -0.8   0.7  -1.3   2.9   1.8  -0.3  -0.1
## 2004   1.7  -0.2  -0.5   1.8   0.0   0.3   1.3  -1.2   0.2  -1.6   1.0  -1.4
## 2005   2.4  -0.2   0.4   0.5  -0.7   0.3  -0.9   0.1   1.0  -0.5   1.7   0.2
## 2006   0.5  -0.1   0.6   2.2  -1.6   1.4   0.4   0.4   0.0   0.9   0.4   0.3
## 2007  -1.3   0.7   0.0  -0.5   1.3   0.1   0.1   3.0  -1.9   2.2  -1.7   0.8
## 2008   0.6   0.3  -1.1  -0.7   0.6  -3.3  -0.2  -4.0   1.3  -2.7  -7.3  -8.4
## 2009  -8.3  -7.3   1.0   3.5   3.0   1.6   1.2   1.3   3.0   2.4   2.0   0.8
## 2010   2.5   0.4   0.2   1.1  -0.2  -0.8   1.1   0.4   1.6  -2.9   1.6   0.6
## 2011   0.5   0.6 -17.2   1.9   6.1   4.0   1.2   1.7  -0.9   1.8  -2.2   2.0
## 2012   0.4  -0.2  -0.2  -0.5  -1.8  -0.8  -0.5  -1.5  -2.1   0.3  -0.9   1.3
## 2013  -1.6   1.7   1.2   0.0   1.6  -1.1   1.6   0.2   1.0   0.2   0.6   0.0
## 2014   2.0  -1.1   1.5  -4.6   2.3  -1.6  -0.2  -0.6   1.2  -0.3   0.0  -0.5
## 2015   3.0  -3.1  -0.5   0.2   0.0   0.9  -0.1  -1.7   2.0   0.1  -0.8  -1.4
## 2016   1.6  -0.9   0.5  -0.4  -0.8   0.7   0.6   0.7   0.2   0.3   1.0   0.0
## 2017  -1.1   0.7  -0.1   2.6  -1.8   1.0  -0.8   1.5  -1.0   0.3   0.9   1.6

# 階差系列の平均
mean_dx <- mean(dx)
mean_dx

## [1] 0.08935282

階差の計算は1978年2月から始まるので、1978年1月の観測値が欠如している。

# 時系列プロット：階差系列
par(family = "HiraginoSans-W3") #日本語が文字化けしないようにおまじない
plot(dx, xlab = '年月', ylab = '鉱工業生産指数の階差系列')
abline(a = mean_dx, b = 0, col = 'blue') #平均を青線で描画

２つの成長率

離散型成長率 \[ r_t = \frac{X_t}{X_{t-1}} - 1 \]
連続型成長率（株式の成長率など観測個数の多い系列に用いる） \[ \Delta \log(X_t) = \log(X_t) - \log(X_{t-1}) \]
成長率が小さな値を取るときは、２つの成長率には違いが見られず、ほとんど同じ値になる。

#成長率の計算
n <- length(x)
r_d <- x[2:n] / x[1:n-1] - 1 #離散型成長率
r_c <- log(x[2:n]) - log(x[1:n-1]) #連続型成長率

# 時系列プロット：成長率
par(family = "HiraginoSans-W3") #日本語が文字化けしないようにおまじない
plot(r_d, type = 'l', xlab = '年月', ylab = '成長率')
lines(r_c, lty = 2) #点線で連続型成長率を重ね書き

２つの成長率にはほとんど差がない。

自己相関

標本自己相関係数

自己相関係数（auto-correlation coefficient：AC）
- 観測値の時間的な依存度を測る尺度
標本AC係数
- 原系列をある期間ずらした系列の間で計算された標本相関係数
- 例：１次の標本AC係数
  - 原系列$\{X_1, X_2, \ldots, X_n\}$と１期ずらした系列（ラグ系列という）$\{\mbox{なし}, X_1, X_2, \ldots, X_{n-1}\}$の間の相関係数
  - 観測個数が１個減ることに注意
- １期ずつ期間をずらして求めた標本AC係数をずれ幅の順に整理したものを標本AC関数と呼ぶ。
ラグ期間が$k$である$k$次の標本AC係数 \[ \begin{align} AC(k) &= \frac{\sum_{t=k+1}^{n} (X_t - \bar{X})(X_{t-k} - \bar{X})}{\sqrt{\sum_{t=1}^{n} (X_t - \bar{X})^2 \sum_{t=k+1}^{n} (X_{t-k} - \bar{X})^2}} \\ &\fallingdotseq \frac{\sum_{t=k+1}^{n} (X_t - \bar{X})(X_{t-k} - \bar{X})}{\sum_{t=1}^{n} (X_t - \bar{X})^2} \end{align} \]
- $\sum_{t=1}^{n} (X_t - \bar{X})^2$と$\sum_{t=k+1}^{n} (X_{t-k} - \bar{X})^2$は大まかには同じ
  - 「$\fallingdotseq$」はほぼ等しい（ニアリーイコール）という意味

母自己相関係数

標本AC係数の背後にあるのは母自己相関係数（真のAC係数）で、理論的な値である母自己相関係数を推定する方法が標本AC係数である。

真の値は次のように定義される。 \[ \rho(k) = \frac{Cov(X_t, X_{t-k})}{V(X_t)} \]
- 分子は母自己共分散、分母は母分散
$X_t$と$X_{t-k}$独立に分布するなら、分子の母自己共分散は０になるから母自己相関係数も０（後述のホワイト・ノイズ）

理論的な自己共分散や分散は、十分に$t$が大きければ、観測時点$t$に依存しないという前提条件が置かれる。

母自己共分散は時間差$k$のみに依存し、分散は定数になる。
このような条件を定常性の条件という。

標本自己相関係数の分散

データ分析では、「$k$次の母自己相関係数が０」であるという帰無仮説を検定することも必要となる。

「$k$次の自己相関が０」という帰無仮説のもとで、$k$次の標本自己相関係数の分散は \[ V(AC(k)) = \frac{1}{n} \left( 1 + 2\sum_{j=1}^{k-1}\rho(j)^2 \right) \] にもとづいて推定できる（$k-1$次までの母自己相関は０とされない）
$\rho(j)$を$AC(j)$に代えてこの値を推定する
１次の母自己相関係数が０という帰無仮説の検定では$V(AC(1))=\frac{1}{n}$となる。
有意水準が5%なら、$AC(1) \pm 2 \sqrt{\left(\frac{1}{n}\right)}$（いわゆる２シグマ）が有意であるかどうかを検定するための境界を示す。

# 鉱工業生産指数の原系列の標本AC関数と2シグマ
acf(x, lag.max = 20, ci.type = "ma")

# 鉱工業生産指数の連続型成長率の標本AC関数と2シグマ
r <- diff(log(x)) #連続型成長率
acf(r, lag.max = 20, ci.type = "ma")

横軸のLagの単位は1年（12ヶ月分）となっていることに注意。
原系列は時点が離れるほど相関は低くなっているが、20次以降もAC係数が２シグマを超えており有意な相関がある。
連続型成長率は2次のAC係数が２シグマを超えており有意であるが、それ以降の相関は無視できる。

自己回帰（AR）法

線形回帰分析では、被説明変数$X_t$の値が経済的な意味を備えた回帰式によって定まる（回帰式が経済構造の骨組みを示す）と考える。

このような考え方は構造的アプローチといわれる。

一方、時系列分析では被説明変数$X_t$の変化は、同じ変数の過去の値（例えば$X_{t-1}$）によって説明されるとする。

説明変数が被説明変数の過去の値（ラグ値）になっている回帰式を、自己回帰式（auto-regression：AR）と呼ぶ。

１次の自己回帰式

簡単な１次の自己回帰式（AR(1)と表記する）は次のようになる。 \[ X_t = c + \phi X_{t-1} + \epsilon_t \]
- この式では、被説明変数$X_t$の変動が、１期前の観測値$X_{t-1}$によって説明されると考えられている。
- $\epsilon_t$は回帰式の誤差項に似ており、$X_{t-1}$よっては説明しえない変動をまとめて示す。
  - 時系列分析ではノイズ（noise）、ショック（shock）、あるいはイノベーション（innovation）と呼ぶ。
- ノイズ$\epsilon_t$は、平均が０、分散は一定、かつ互いに独立に分布する確率変数であると仮定される。
  - このような標準的なノイズをホワイト・ノイズと呼ぶ。
- 変数の過去の値$X_{t-1}$をラグつき変数という。
- 回帰係数$\phi$を１次の自己回帰係数と呼ぶ。
  - $|\phi| < 1$という（定常性の）条件が課せられる。
経済分析において自己回帰式は、分析の対象となっている変数の変動を近似するために利用される。
- 例えば、GDPは、さまざまな要因が未知で複雑なメカニズムに投入投入（インプット）されて産出（アウトプット）されたものと理解される。
  - メカニズムを簡潔な式の体系で表現しようとするのが構造式アプローチである。
  - 時系列分析では、このメカニズムをあまりに複雑で捉えることができないブラック・ボックスとみなし、アウトプットだけを捉え、その変動を近似しようと試みる。
  - このような単純化のおかげで経済変数に関する時系列予測は利用可能になる一方で、経済構造を無視するという批判を受けることもある。

AR過程のノイズによる表現

１次の自己回帰式は、定数項を分解して次のように書くこともできる。 \[ X_t - \mu = \phi (X_{t-1} - \mu) + \epsilon_t \]
- 定数項は$c = (1-\phi)\mu$と定義され、系列$X$の平均$\mu$と自己回帰係数$\phi$で決まる。
- $X_t$の平均と$X_{t-1}$の平均は、定常性により値が変わらないとし、$\mu = \frac{c}{1-\phi}$となる。
$X_t = c + \phi X_{t-1} + \epsilon_t$はすべての$t$について成立するから、$X_0$を０期の値として、$t=1$について \[ X_1 = c + \phi X_0 + \epsilon_1 \]
- $X_0$を初期値という。
第２期の変数は \[ X_2 = c + \phi X_1 + \epsilon_2 = (1+\phi)c + (\epsilon_2 + \phi \epsilon_1) + \phi^2 X_0 \]
$|\phi| < 1$とし、同じ操作を繰り返していくと、第$t$期は \[ X_t = \mu + (\epsilon_t + \phi \epsilon_{t-1} + \cdots + \phi^{t-1}\epsilon_1) + \phi^t X_0 \]
- $\mu$は$X_t$の平均で（$t$が十分大きければ） \[ \mu = (1 + \phi + \cdots + \phi^{t-1})c \fallingdotseq \frac{c}{1-\phi} \]

phi <- 0.9
t <- 0:20
# \phiのt乗の推移(t=0,1,...,20)
plot(t, phi^t)

$t$が十分大きければ$\phi^t$を０とみなしてよいから、$X_t$は次のように表現できる。 \[ \begin{align} X_t &= \mu + u_t \\ u_t &= \epsilon_t + \phi \epsilon_{t-1} + \phi^2 \epsilon_{t-2} + \cdots \end{align} \]
- $X_t$は平均$\mu$を除いて、無限次のホワイト・ノイズの系列によって生成されていることがわかる。
この式を自己回帰過程の移動平均化と呼ぶ。
- 変数がノイズのみで表されているので、エラーショック表現ともいわれる。
- １次より大きい一般の自己回帰式も、定常性の条件が満たされていれば、このような無限次移動平均あるいはエラーショック表現を導くことができる。
  - AR(1)では$|\phi|<1$が定常性の条件になる。
  - 定常性の条件は、時系列変数の理論的な平均や分散が観測時点に依存しないための条件であり、また他方で、自己回帰式を無限次の移動平均で表現するための条件にもなっている。

高次の自己回帰式

$p$次式AR($p$)は \[ \begin{align} X_t &= c + \phi_1 X_{t-1} + \phi_2 X_{t-2} + \cdots + \phi_{p}X_{t-p} + \epsilon_t \\ &\quad (t = p+1, p+2, \ldots, n) \end{align} \]

この式は最小２乗法で推定できるが、後述のMAやARIMAと共に最尤法で推定されることも多い。

推定における基本的な注意

初期値に観測値が$p$個必要で、観測個数は$p$個減少する。
式のフィット（適合度）のよさは、各自己回帰係数に関する$t$値と、式全体の決定係数$R^2$を用いて検討する。
$t$検定により、特定の項$X_{t-s}$を排除する場合、$s$次より低次のラグはすべて残し、高次のラグはすべて除くように式を定める。 \[ \begin{align} X_t &= c + \phi_1 X_{t-1} + \phi_2 X_{t-2} + \cdots + \phi_{p}X_{t-(s-1)} + \epsilon_t \\ &\quad (t = s, s+1, \ldots, n) \end{align} \]
- $X_{t-s}$は除去するが高次の項である$X_{t-(s+1)}$は残す、といった判断はしないのが普通
残差に関する検定では、次節で説明される残差に関する自己相関関数、ふろしき検定などを使う。

推定例と次数の選択

鉱工業生産指数の連続型成長率について自己回帰式を推定する。

まず、4次の自己回帰式AR(4)を推定する。

# AR(4)の推定
ar4 <- arima(r, order = c(4, 0, 0))
ar4

## 
## Call:
## arima(x = r, order = c(4, 0, 0))
## 
## Coefficients:
##           ar1     ar2     ar3      ar4  intercept
##       -0.0129  0.1343  0.0514  -0.0006     0.0011
## s.e.   0.0457  0.0457  0.0456   0.0457     0.0010
## 
## sigma^2 estimated as 0.0003234:  log likelihood = 1245.09,  aic = -2478.19

Rではarima関数を使って簡単にARモデルを最尤法により推定できる。
- orderの設定で後述のMAモデルやARIMAモデルも推定できる。
- オプションでmethod = "CSS"とすれば最小二乗法（正確には条件付き最小二乗法）による推定もできる。
Coefficeintsには各ラグ項の係数とともに、係数の標準誤差s.e.が表示される。
誤差項$\epsilon_t$の分散sigma^2の推定値は0.0003234で、対数尤度log likelihoodは1245
aicは赤池情報量基準（Akaike information criterion：AIC）」と呼ばれるモデルの総括的な比較基準で、最小値を与えるモデルが最も望ましいと考えられる。 \[ AIC = -2 \times \mbox{対数尤度} + 2 \times \mbox{パラメータ数} \]
- パラメータ数は$\mbox{ラグ項}+\mbox{定数項}+\mbox{誤差項の分散}$として計算する
  - AR(4)では$4+1+1=6$
係数の推定値を標準偏差で割ったもの（の絶対値）が$t$値で、各係数が0とう帰無仮説を検定するために使われる。
- 線形回帰の場合と異なり、帰無仮説のもとでの$t$統計量の分布は、近似的に標準正規分布になる。
- 簡単化のために$t$値が2より大きければ有意と考えることが多い。（理由は？考えてみよう）
- 2次の係数（ar2）以外は有意にならないことがわかる。

高次の項が有意でないので、字数を下げて推定する。

# AR(3)の推定
ar3 <- arima(r, order = c(3, 0, 0))
ar3

## 
## Call:
## arima(x = r, order = c(3, 0, 0))
## 
## Coefficients:
##           ar1     ar2     ar3  intercept
##       -0.0129  0.1342  0.0514     0.0011
## s.e.   0.0456  0.0452  0.0456     0.0010
## 
## sigma^2 estimated as 0.0003234:  log likelihood = 1245.09,  aic = -2480.19

# AR(2)の推定
ar2 <- arima(r, order = c(2, 0, 0))
ar2

## 
## Call:
## arima(x = r, order = c(2, 0, 0))
## 
## Coefficients:
##           ar1     ar2  intercept
##       -0.0060  0.1340     0.0011
## s.e.   0.0453  0.0453     0.0009
## 
## sigma^2 estimated as 0.0003242:  log likelihood = 1244.46,  aic = -2480.92

# AR(1)の推定
ar1 <- arima(r, order = c(1, 0, 0))
ar1

## 
## Call:
## arima(x = r, order = c(1, 0, 0))
## 
## Coefficients:
##           ar1  intercept
##       -0.0067     0.0011
## s.e.   0.0457     0.0008
## 
## sigma^2 estimated as 0.0003302:  log likelihood = 1240.13,  aic = -2474.26

AICを比較する。

AIC(ar1, ar2, ar3, ar4)

##     df       AIC
## ar1  3 -2474.256
## ar2  4 -2480.918
## ar3  5 -2480.188
## ar4  6 -2478.188

最小AIC基準によればAR(2)が選ばれる。

自己回帰の次数を決める際には、各式の推定結果だけでなく、残差の性質を検討する必要がある。

特に残差に系列相関が残っていない式を選択することが必要であり、そのために様々な検討を加える。
ここでは、系列r、AR(1)によりフィットした値、rとフィット値との差である残差の時系列プロットを見る。

# 成長率とフィット値のプロット
par(family = "HiraginoSans-W3") #日本語が文字化けしないようにおまじない
ts.plot(r, r - ar2$residuals, gpars=list(xlab="年月", ylab="成長率", lty=c(1:2), col=c(1:2))) #成長率とフィット値

# 残差のプロット
par(family = "HiraginoSans-W3") #日本語が文字化けしないようにおまじない
ts.plot(ar2$residuals, gpars=list(xlab="年月", ylab="残差", lty=3, col=3))

2009年頃の世界金融危機と2011年の大震災の影響により、残差が大きな値を示している。
- 誤差項$\epsilon_t$の分散sigma^2の推定値は0.0003302で、２シグマ0.036となる。
- 上記の2つの時期で、残差は２シグマを大きく超えている。
より詳しい検討方法は後述。

標本偏自己相関（PAC）係数

自己回帰式の推定において、最高次の項の係数推定値は、標本偏自己相関（partial auto-correlation：PAC）係数と呼ばれる。

自己回帰式の望ましい次数を選ぶ際に、PACが重要な役割を果たす。
個々のラグ変数の$t$値よりも、**最高次係数の$t$値$をその判断材料にする。
- 最高次のラグが有意であれば、それ以下のラグ次数が有意でなくとも、AR式に低次の項を残しておく。

回帰式は説明変数が少なすぎると推定に偏り（バイアス）を生じる。

説明変数が多すぎると、回帰式に推定上の偏りはないが、推定に無駄が生じ、有効性が失われる。
観測期間（サンプルサイズ）が大きければ、偏りを避けるために少々次数を高くしておいた方が安心。

PACの計算法

$p$次のPAC、PAC($p$)は、AR($p$)式の$p$次のラグ項係数の最小２乗推定値として求めることができる。

Rでは、関数pacfを用いて計算することができる。

# PAC
pacf(r, lag.max = 20)

２シグマ線（点線）は、すべての自己相関が０であるという帰無仮説の元で計算されている。
後述の診断検定に伴う残差のPACでは、PACの分散は$1/n$であり、線形回帰からもとまる係数の標準偏差と異なる値になる。
- $n=480-1=479$とすると、標準偏差は0.0457となり、AR(1)の１次の係数の標準偏差0.0457と（たまたま）同じになっている。
図の２シグマの値は0.0914であり、3次以降は有意でない。
- ここでは、PACの観点からも、AR(2)が望ましい式と考えられる。

参考文献

計量経済学（有斐閣）http://www.yuhikaku.co.jp/books/detail/9784641053854
基礎からわかる時系列分析（技術評論社）https://gihyo.jp/book/2018/978-4-7741-9646-6
経済・ファイナンスデータの計量時系列分析（朝倉書店）https://www.asakura.co.jp/books/isbn/978-4-254-12792-8/

経営のための統計学II

第3回時系列分析1

高橋慎

2021/10/08

時系列データの種類

時系列プロット

鉱工業生産指数データ（csvファイル）の読み込み

Rにおける時系列データの扱い

原系列

対数系列

階差系列

２つの成長率

自己相関

標本自己相関係数

母自己相関係数

標本自己相関係数の分散

自己回帰（AR）法

１次の自己回帰式

AR過程のノイズによる表現

高次の自己回帰式

推定例と次数の選択

標本偏自己相関（PAC）係数

PACの計算法

参考文献

経営のための統計学II

第3回 時系列分析1

高橋 慎

2021/10/08

時系列データの種類

時系列プロット

鉱工業生産指数データ（csvファイル）の読み込み

Rにおける時系列データの扱い

原系列

対数系列

階差系列

２つの成長率

自己相関

標本自己相関係数

母自己相関係数

標本自己相関係数の分散

自己回帰（AR）法

１次の自己回帰式

AR過程のノイズによる表現

高次の自己回帰式

推定例と次数の選択

標本偏自己相関（PAC）係数

PACの計算法

参考文献

第3回時系列分析1

高橋慎