「続・高校生の統計学」のプログラムコードを掲載しています。

RstudioやPosit Cloud.
GoogleColabのR環境.

においてコピー＆ペーストで利用できます。
初めは指導者のもとで行うのが良いでしょう。
（電子書籍はデジタル著作権管理（DRM）仕様のためコピーできません）

インストールせずにブラウザでプログラムを実行するには、次の方法があります。

Posit Cloud で Rを利用（無料版があります）.
参考サイト Rstudio_Cloudの利用.
GoogleColab でランタイムのタイプを R に変更して利用.
参考サイト Google-Colabの利用.

以下、電子書籍との章番号を対応させるため、項目名は残してあります。

必要なパッケージを用いるのに必要なコードで書籍に表示されていないものもありますが、以下では表示してあります。
例えば、library(mvtnorm)、library(scatterplot3d).

1 はじめに.

2 カイ二乗検定

2.1 四分表

2.1.1 独立.

データの出力には、多項乱数
　　　rmultinorm(乱数の数, 調査数, 比率)　　を用います。

＜注意＞.
簡単のため、行の周辺度数を固定して考えます。
各行の中での効果の人数は二項分布に従い、考えやすくなります。
データ数が十分大きければ二項分布は正規分布で近似できるため、探究的学習の可能性が広がります。

N=500、$\pi_1=3/5$ とすると、行の比率は $3:2$ なので、
　　　１行目の合計は300、２行目の合計は200　となります。

１行目の合計 300 を、$\pi_2=3/10$ に従って $3:7$ に分けると、

ab <- rmultinom(n=1, size = 300, prob = c(3/10,7/10)) 
t(ab)#転置

##      [,1] [,2]
## [1,]   84  216

新薬を300人に投与したとき、(効果のあった人数, なかった人数) です。
１番目の数字が a 、２番目の数字が b になります。

２行目の合計 200 を、$\pi_2=3/10$ に従って $3:7$ に分けると、

cd <- rmultinom(n=1, size = 200, prob = c(3/10,7/10)) 
t(cd)#転置

##      [,1] [,2]
## [1,]   64  136

ブラセポ（偽薬）を200人に投与したとき、(効果のあった人数, なかった人数) です。
１番目の数字が c 、２番目の数字が d になります。

四分表の形式で表示するには、matrix()関数を使います。

a <- ab[1,];b <- ab[2,];c <- cd[1,];d <- cd[2,]
t(matrix(c(a,b,c,d),2,2))#転置

##      [,1] [,2]
## [1,]   84  216
## [2,]   64  136

これで ad-bc が1つ計算されます。
n=1000 とすれば、1000個の ad-bc が計算され、ヒストグラムを描くことができます。

ab <- rmultinom(n=1000, size = 300, prob = c(3/10,7/10)) 
cd <- rmultinom(n=1000, size = 200, prob = c(3/10,7/10)) 
a <- ab[1,];b <- ab[2,];c <- cd[1,];d <- cd[2,]
x <- a*d-b*c
hist(x,prob=T,col="lightblue")

正規分布に従っているかどうかをみるため、データをscale()関数で標準化したヒストグラムと標準正規分布の理論曲線を比べます。

hist(scale(x),prob=T,col="lightblue")
curve(dnorm(x, 0,1),add=T,,col="blue")

統計量 $X=ad-bc$ は正規分布に従っているようです。

2.1.2 ad-bc の分散

独立という仮説のもとで $X=ad-bc$ は平均０の正規分布に従うとして、
$X$ の分散が分かれば、標準化したものを統計量とすることができます。

ピアソンの統計量は、

　　　$T=\frac{(ad-bc)^2N}{(a+b)(c+d)(a+c)(b+d)} \sim \chi^2(1)$

$\chi^2(1)$ は標準正規分布の２乗なので、

　　　$t=\frac{(ad-bc)\sqrt{N}}{\sqrt{(a+b)(c+d)(a+c)(b+d)}} \sim N(0,1)$

$X=ad-bc$ の分散は、
　　　$v=\frac{(a+b)(c+d)(a+c)(b+d)}{N}$ 　　となりそうですが、
このシミュレーションモデルでは $v$ は変動する量です。

行の周辺度数 $(a+b)=300,~ (c+d)=200$ は定数ですが、列の周辺度数は変動するからです。

代替として、列の周辺度数を、 $500\times \frac{3}{10}=150,~~500\times \frac{7}{10}=350$　と考え、
　　　$V=\frac{300\times 200 \times 150\times 350}{500}$ 　　としてみましょう。
　　　

m <- c(mean(a+c),mean(b+d)) #列の周辺度数
m

## [1] 149.933 350.067

v <- c(300*200*150*350/500,mean(300*200*(a+c)*(b+d)/500),var(x))
v

## [1] 6300000 6286513 6455414

列の周辺度数は、ほぼ指定した比率通りです。また、代替値で求めた分散は実際の分散に近くなっています。

2.1.3 ピアソンの予想

筆者の想像ですが、ピアソンは上記の例で
$X=ad-bc$ の分散を　$V=\frac{300\times 200 \times 150\times 350}{500}$　　と考え、
　　　$z=\frac{ad-bc}{\sqrt{V}}\sim N(0,~1)$　と予想したのではないでしょうか。
実際、

V <- 300*200*150*350/500
z <- x/sqrt(V)
hist(z,prob=T,col="lightblue")
curve(dnorm(x, 0,1),add=T,col="blue")

数理的証明は難しそうですが、標準正規分布に近いようです。

2.1.4 数理的考察

2.2 独立検定

対立仮説「新薬に効果がある」のもとでシミュレーションを行い、独立性を検定できるか確かめます。

2.2.1 z検定

帰無仮説のもとでの理論曲線を青で表示しておきます。

N=500
ab <- rmultinom(n=1000, size = 300, prob = c(3/10,7/10)) 
cd <- rmultinom(n=1000, size = 200, prob = c(2/10,8/10)) 
a <- ab[1,];b <- ab[2,];c <- cd[1,];d <- cd[2,]
x <- a*d-b*c
t <- x*sqrt(N)/sqrt((a+b)*(c+d))/sqrt((a+c)*(b+d))
hist(t,prob=T)
curve(dnorm(x, 0,1),add=T,col="blue")
abline(v = qnorm(0.95,0,1),col="red")

赤の縦線は帰無仮説のもとでの棄却境界であり、ヒストグラムで赤の縦線の右側の値が得られたとき帰無仮説は棄却されます。つまり２群の比率の差が検出されます。
この例では調査数が大きいので、わずかの比率の差でも検出できるようです。

2.2.2 カイ二乗検定

Z値の二乗であるカイ二乗値を用いて、比率の差を検定します。

N=500
ab <- rmultinom(n=1000, size = 300, prob = c(3/10,7/10)) 
cd <- rmultinom(n=1000, size = 200, prob = c(2/10,8/10)) 
a <- ab[1,];b <- ab[2,];c <- cd[1,];d <- cd[2,]
x <- a*d-b*c
k <- x^2*N/(a+b)/(c+d)/(a+c)/(b+d)
hist(k,prob=T)
curve(dchisq(x,1),add=T,col="blue")
abline(v = qchisq(0.95,1),col="red")

ヒストグラムで赤の縦線の右側の値が得られたとき帰無仮説は棄却されます。

2.3 適合度検定

適合度検定とは「分割表の観測度数が理論比率にあっているか」を検証する検定方法です。
検定統計量は　$\displaystyle \chi_0^2=\sum_{i=1}^{k}\frac{(o_i-e_i)^2}{e_i}$（理論比率からの乖離度を表す統計量）

ピアソンは、どのようにしてこの検定量にたどり着いたのでしょうか？
四分表の独立性を評価する検定量と、理論度数からの乖離度を表す統計量を比較することで調べてみましょう。

2.3.1 乖離度(適合度)

2.3.2 ピアソンの定理

2.3.3 二項分布.

二項分布 $X \sim B(10,0.5)$に従う二項乱数を100個出力し、データと理論分布を比較します。

a <- rbinom(100,10,0.5)#二項乱数
h <- hist(a,breaks=seq(-0.5, 10.5, by = 1),plot = FALSE)
par(family = "HiraKakuProN-W3") #日本語フォントの指定
par(mfrow=c(1,2))
barplot(h$counts,names.arg=h$mids,ylim=c(0,30),main="データ分布",col="lightblue")
x=0:10
m=dbinom(x,10,0.5)*100
barplot(m,names.arg=c(0:10),ylim=c(0,30),main="理論分布",col="lightblue")

適合度検定量 $\chi_0^2$ を計算します。

x <- h$counts
K <- sum((x-m)^2/m)
K

## [1] 9.990603

ピアソンの定理によれば、
$N=100$、各区分の確率を　$p_0,~p_1,~\cdots p_{10}$ とすると、
n区分の適合度検定量は自由度(n-1)のカイ二乗分布に従います。即ち、

　　　$T=\displaystyle \sum_{i=0}^{10}\frac{(X_i-Np_i)^2}{Np_i}\sim \chi^2(10)$　

これをシミュレーションで確かめましょう。

2.3.3.1 シミュレーション.

1000回繰り返して、自由度１０のカイ二乗分布の理論曲線と比較します。

K <- c()
for(i in 1:1000){
  a <- rbinom(100,10,0.5)#二項乱数
  h <- hist(a,breaks=seq(-0.5, 10.5, by = 1), plot = FALSE)
  x <- h$counts
  m <- dbinom(0:10,10,0.5)*100
  k <- sum((x-m)^2/m)
  K <- c(K,k)
}
hist(K,xlim=c(0,20),prob=T,breaks="scott",col="lightblue")
curve(dchisq(x,10),xlim=c(0,20),add=T,col="blue")

mean(K)

## [1] 9.687438

mean(K) の値から見て、自由度１０のカイ二乗分布に従いますが、「理論分布と少しずれがあるように感じられます。

2.3.3.2 コクランの規則.

2.3.3.3 補正シミュレーション.

右端の３区分と、左端の３区分を合併し、全部で７区分とすれば、適合度検定の自由度は 7-1=6 となるはずです。

K <- c()
for(i in 1:1000){
  a <- rbinom(100,10,0.5)#二項乱数
  h <- hist(a,breaks=seq(-0.5, 10.5, by = 1), plot = FALSE)
  x <- c(sum(h$ counts[1:3]),h$ counts[4:8],sum(h$ counts[9:11]))
  m <- dbinom(0:10,10,0.5)
  n <- c(sum(m[1:3]),m[4:8],sum(m[9:11]))*100
  k <- sum((x-n)^2/n)
  K <- c(K,k)
}
hist(K,xlim=c(0,20),prob=T,breaks="scott",col="lightblue")
curve(dchisq(x,6),xlim=c(0,20),add=T,col="blue")

mean(K)

## [1] 5.936908

かなり近似が良くなっており、自由度は６です。

2.3.3.4 推定母数があるとき.

期待度数を求めるために観測値からr個の母数を推定していれば、自由度は r 減ります。

　　　$T=\displaystyle \sum_{i=1}^{k}\frac{(X_i-Np_i)^2}{Np_i}\rightarrow \chi^2(k-r-1)$

p=0.5 の代わりに、推定値 p = mean(a)/10 を用います

理論累積確率 dbinom(h$breaks, 10, 0.5)　のところを.
推定値を用いて　dbinom(h$breaks, 10, mean(a)/10).

としてシミュレーションすると自由度は１減り、自由度５となるはずです。

K <- c()
for(i in 1:1000){
  a <- rbinom(100,10,0.5)#二項乱数
  h <- hist(a,breaks=seq(-0.5, 10.5, by = 1), plot = FALSE)
  x <- c(sum(h$ counts[1:3]),h$ counts[4:8],sum(h$ counts[9:11]))
  m <- dbinom(0:10,10,mean(a)/10)
  n <- c(sum(m[1:3]),m[4:8],sum(m[9:11]))*100
  k <- sum((x-n)^2/n)
  K <- c(K,k)
}
hist(K,xlim=c(0,20),prob=T,col="lightblue")
curve(dchisq(x,5),xlim=c(0,20),add=T,col="blue")

mean(K)

## [1] 4.901634

統計量の平均は５となっており、自由度が５となっていることが分かります。
自由度５のカイ二乗分布の理論曲線はヒストグラムによく当てはまっています。

3 時系列分析.

3.1 AR(1)モデル.

例　$y_1=3,~~~y_{t}=3-0.5 y_{t-1}+\varepsilon _{t}$

yt=3 #初期値
Yt=c(3)
for (i in 1:200) {
  yt <- 3-0.5*yt+rnorm(1,0,1)
  Yt <- c(Yt,yt)
}
plot(Yt,type="l",xlab="time")
abline(h=mean(Yt),col="red")
abline(h=mean(Yt)+1.96*sd(Yt),col="blue", lty=2)
abline(h=mean(Yt)-1.96*sd(Yt),col="blue", lty=2)

上記の例では、

$\mu=\frac{c}{1-\phi}=\frac{3}{1+0.5}=2$
$V=\frac{\sigma^2}{1-\phi^2}=\frac{1}{1-(-0.5)^2}\fallingdotseq1.33~~~\therefore sd\fallingdotseq1.15$

mv <- c(mean(Yt),sd(Yt))
mv #平均と標準偏差

## [1] 1.970103 1.086234

実データの平均と標準偏差は理論値とほぼ一致しています。

3.1.1 自己相関係数.

自己相関係数は、過去データが現在に及ぼす影響の度合いを表します。

１過去（１つ前）データとの自己相関係数 $\rho_1$ は、
　　　$\rho_1=Cor(Y_t,~Y_{t-1})$　ですが、
　　　$\rho_1=Cor(Y_t,~Y_{t+1})$　でも計算することができます。

上記の例で、 $Y_t$ から最初の１項を除けば、$Y_{t+1}$ になります。
最初の１項を除く操作は、tail(Yt,-1) です。はじめの６項を head()関数で表示すると、

head(Yt)

## [1] 3.0000000 0.4483732 2.4536351 0.8053067 2.4708649 1.0340166

head(tail(Yt,-1))

## [1] 0.4483732 2.4536351 0.8053067 2.4708649 1.0340166 3.2922107

tail(Yt,-1)によって、項が１つスライドしていることが分かります。

データの個数を揃えるため head(Yt,-1) によってYt の最後の１項を除き、相関係数を計算します。

r1 <- cor(head(Yt,-1),tail(Yt,-1))
r1

## [1] -0.4884666

これは、ラグ１（ズレ１）の標本自己相関係数で、漸化式の係数に近い値となっています。

3.1.2 コレログラム.

Rにはacf関数(autocorrelation function) が用意されており、自己相関係数を計算し図示することができます。

plot=F で数値を表示し、
plot=T でグラフを表示します.

グラフをコレログラムと呼びます。

#自己相関係数
acf(Yt,plot=F,lag.max=7)

## 
## Autocorrelations of series 'Yt', by lag
## 
##      0      1      2      3      4      5      6      7 
##  1.000 -0.487  0.162  0.020 -0.059  0.081 -0.046 -0.072

acf(Yt,plot=T)

青の点線で挟まれた範囲は、無相関という帰無仮説が棄却されない範囲です。この帯の間にグラフがあれば、そのラグで相関はないと判断されます。
上図では、はっきりした影響を及ぼすのは１過去までですが、２過去も若干影響が残っているようです。

acf(Yt,plot=T,type="p")

間接的影響を除くと、直接影響を与えるのは１過去のデータのみであることが分かります。

3.1.3 漸化式との関係.

誤差項のない隣接２項間の漸化式は、
　　　$y_{t}=3-0.5 y_{t-1}$
一般項は漸化式を解くことにより、
　　　$y_{t}=2+(y_0-2)(-0.5)^t$
変化の様子は、

t=1:10
y=2+(3-2)*(-0.5)^t
y

##  [1] 1.500000 2.250000 1.875000 2.062500 1.968750 2.015625 1.992188 2.003906
##  [9] 1.998047 2.000977

plot(y,type="l",ylim=c(1.8,2.2),xlab="time")
abline(h=2,col="lightblue")

漸化式の一般項は２に収束し、変動もシンプルです。
時系列では誤差項が加わるため、変動が複雑化していることが分かります。

3.2 定常性.

3.2.1 定常 AR(1).

3.2.1.1 $c=2,~~\phi=0.5$ のとき.

　$y_1=3,~~~y_{t}=2+0.5 y_{t-1}+\varepsilon _{t}$.

yt=3
Yt=c(3)
for (i in 1:(200-1)) {
  yt <- 2+0.5*yt+rnorm(1,0,1)
  Yt <- c(Yt,yt)
}
head(Yt)

## [1] 3.000000 4.453120 3.318068 2.591920 4.101439 1.870887

length(Yt)

## [1] 200

plot(Yt,type="l",xlab="time")
abline(h=4,col="red")
abline(h=mean(Yt)+1.96*sd(Yt),col="blue", lty=2)
abline(h=mean(Yt)-1.96*sd(Yt),col="blue", lty=2)

＜平均と標準偏差＞.
　実データによる標本値と理論値を比較します。

mv <- c(mean(Yt),sd(Yt))
mv #標本値

## [1] 4.045303 1.197966

MV <- c(2/(1-0.5),sqrt(1/(1-0.5^2)))  
MV #理論値

## [1] 4.000000 1.154701

標本値と理論値はほぼ一致しています。

3.2.2 定常の意味.

3.2.2.1 $c=2,~~\phi=0.5$ のとき.

　$y_1=3,~~~y_{t}=2+0.5 y_{t-1}+\varepsilon _{t}$

yt=3
Yt=c(3)
for (i in 1:(100-1)) {
  yt <- 2+0.5*yt+rnorm(1,0,1)
  Yt <- c(Yt,yt)
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(Yt,type="l",xlab="time",main="誤差項あり")
abline(h=4,col="red")

yt=3
Yt=c(3)
for (i in 1:(100-1)) {
  yt <- 2+0.5*yt
  Yt <- c(Yt,yt)
}
plot(Yt,type="l",xlab="time",main="誤差項なし")
abline(h=4,col="red")

時系列データは単調かつ急速に収束します。
誤差項を加えることで、時系列データは平均を挟んだ変動となっています。

3.2.2.2 $c=2,~~\phi=-0.5$ のとき.

　$y_1=3,~~~y_{t}=2-0.5 y_{t-1}+\varepsilon _{t}$

yt=3
Yt=c(3)
for (i in 1:(100-1)) {
  yt <- 2-0.5*yt+rnorm(1,0,1)
  Yt <- c(Yt,yt)
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(Yt,type="l",xlab="time",main="誤差項あり")
abline(h=2/1.5,col="red")

yt=3
Yt=c(3)
for (i in 1:(100-1)) {
  yt <- 2-0.5*yt
  Yt <- c(Yt,yt)
}
plot(Yt,type="l",xlab="time",main="誤差項なし")
abline(h=2/1.5,col="red")

係数がマイナスだと振動しますが、急速に収束します。
誤差項を加えた時系列データは平均を挟んで激しく振動しています。

3.2.3 非定常 AR(1).

3.2.3.1 $c=2,~~~\phi=1.5$ のとき.

　$y_1=3,~~~y_{t}=2+1.5 y_{t-1}+\varepsilon _{t}$
　

yt=3
Y=c(2)
for (i in 1:20) {
  yt <- 2+1.5*yt+rnorm(1,0,1)
  Y <- c(Y,yt)
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(Y,type="l",xlab="time",main="誤差項あり")

yt=3
Y=c(2)
for (i in 1:20) {
  yt <- 2+1.5*yt
  Y <- c(Y,yt)
}
plot(Y,type="l",xlab="time",main="誤差項なし")

誤差項のある時系列も誤差項のない漸化式もデータは発散しています。

3.2.3.2 $c=0,~~~\phi=1$ のとき.

　　$y_{t}=y_{t-1}+\varepsilon _{t}$　、誤差変動：$\varepsilon _{t}\sim N(0,\sigma^2)$　のとき.

　　$y_{t}=\varepsilon _{1}+\varepsilon _{2}+\cdots +\varepsilon _{t}$　（誤差項の和）

これは、ランダムウォークと呼ばれます。時系列も平均値も収束しません。
時系列データと平均変動を並置します。

yt=0
Y=c(0)
M=c()
for (i in 1:50) {
  yt <- 1*yt+rnorm(1,0,1)
  Y <- c(Y,yt)
  M <- c(M,mean(Y))
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(Y,type="l",xlab="time",main="時系列データ")

plot(M,type="l",xlab="time",main="平均変動")

平均もランダムに変動し、一定した平均値を持ちません。

3.2.3.3 $c=2,~~~\phi=1$ のとき.

　　$y_{t}=2+y_{t-1}+\varepsilon _{t}$　、誤差変動：$\varepsilon _{t}\sim N(0,\sigma^2)$　のとき.

　　$y_{t}=2t+\varepsilon _{1}+\varepsilon _{2}+\cdots +\varepsilon _{t}$　（誤差項の和）

この場合、増加トレンドになります。

yt=2
Y=c(2)
M=c()
for (i in 1:50) {
  yt <- 2+1*yt+rnorm(1,0,1)
  Y <- c(Y,yt)
  M <- c(M,mean(Y))
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(Y,type="l",xlab="time",main="時系列データ")
plot(M,type="l",xlab="time",main="平均変動")

3.2.4 まとめ.

3.3 ユール&ウォーカー.

参考サイト.

3.3.1 時系列のモデリング.

3.3.2 ユール・ウォーカー方程式

3.4 MA モデル.

3.4.1 MA(0).

最もシンプルな定常モデルは誤差項だけのモデルです。

　　　$y_{t}=c+\varepsilon _{t}$　　誤差変動：$\varepsilon _{t}\sim N(0,\sigma^2)$

これを M(0) モデルと呼ぶことにします。

c <- 3
y <- rep(0, 200)
e <- rnorm(200)
for (i in 1:length(y)) {
  y[i] <- c + e[i]
}
par(mfrow=c(1,2))
plot(y,type="l",xlab="time",main="Time series")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)
acf(y,plot=T,main="ACF")

mv <- c(mean(y),sd(y))
mv #平均と標準偏差

## [1] 3.015632 1.008150

コレログラムを見ると、M(0)モデルは過去から影響を受けないモデルということができます。
コレログラムに注目すると、変動の特徴が見えてきます。

3.4.2 MA(1).

時刻 $t,~t-1$ でのデータ　$c+\varepsilon _{t},~~c+\varepsilon _{t-1}$ の平均は、

　　　$y_{t}=c+0.5\varepsilon _{t}+0.5\varepsilon _{t-1}$

このモデルによりデータを200個生成します。

c <- 3
y <- rep(0, 200)
y[1] <- 3
e <- rnorm(200)
for (i in 2:length(y)) {
  y[i] <- c + 0.5*e[i]+ 0.5*e[i-1]
}
par(mfrow=c(1,2))
plot(y,type="l",xlab="time")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)
acf(y,plot=T,main="")

mv <- c(mean(y),sd(y))
mv #平均と標準偏差

## [1] 3.0110475 0.6768137

時系列変動だけでは M(0)モデルと区別がつきませんが、コレログラムを見ると M(1)モデルは1過去から影響を受けていることが分かります。
この移動平均モデルは、過去から半分影響を受けるモデルです。

過去の影響のほうが大きいこともあります。加重平均を用いると、様々なモデルを表現できます。
（定数項は変わりません）

　　　$y_{t}=c+0.2\varepsilon _{t}+0.8\varepsilon _{t-1}$
　　　

c <- 3
y <- rep(0, 200)
y[1] <- 3
e <- rnorm(200)
for (i in 2:length(y)) {
  y[i] <- c + 0.8*e[i]+ 0.2*e[i-1]
}
par(mfrow=c(1,2))
plot(y,type="l",xlab="time")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)
acf(y,plot=T)

mv <- c(mean(y),sd(y))
mv #平均と標準偏差

## [1] 2.9677450 0.9423864

コレログラムを見ると、
M(1)モデルは１過去からの影響を受け、それ以前からは影響を受けないモデルということができます。

3.4.3 定常性.

コレログラムを見れば、ARモデルかMAモデルか判断できることになります。

yt=3
Yt=c(3)
for (i in 1:(200-1)) {
  yt <- 2+0.5*yt+rnorm(1,0,1)
  Yt <- c(Yt,yt)
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(Yt,type="l",xlab="time",main="AR(1)モデル")
abline(h=2/0.5,col="red")

y <- rep(0, 200)
y[1] <- 3
e <- rnorm(200)
for (i in 2:length(y)) {
  y[i] <- 3 + e[i]+ 0.5*e[i-1]
}
plot(y,type="l",xlab="time",main="MA(1)モデル")
abline(h=mean(y),col="red")

par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") #日本語フォントの指定
acf(Yt,plot=T,main="AR(1)のコレログラム")
acf(y,plot=T,main="MA(1)のコレログラム")

3.4.4 MA (2).

3.4.4.1 MA(2)の例.

　　　$y_{t}=3+\varepsilon _{t}+5\varepsilon _{t-1}+3\varepsilon _{t-2}$

y <- rep(0, 200)
e <- rnorm(200)
for (i in 3:length(y)) {
  y[i] <- 3 + e[i]  + 5 * e[i - 1] +  3 * e[i - 2]
}
par(family = "HiraKakuProN-W3") #日本語フォントの指定
par(mfrow=c(1,2))
plot(y,type="l",xlab="time",main="MA(2)モデル")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)
acf(y,plot=T,main="コレログラム")

mv <- c(mean(y),sd(y))
mv #平均と標準偏差

## [1] 2.619057 5.948619

3.5 AR(2) モデル.

例　$y_1=0,~y_2=1,~~~y_{t}=2+1.2 y_{t-1}-0.8y_{t-2}+\varepsilon _{t}$

AR(1)では、定常であるためには係数の絶対値が１より小さいことが必要でした。
ここでは、$\phi_1=1.2>1$ ですが、定常性はどうなるでしょうか？

y <- rep(0, 200)
y[1] <- 0
y[2] <- 1
for (i in 3:length(y)) {
  y[i] <- 2 + 1.2 * y[i - 1] -  0.8 * y[i - 2]+rnorm(1,0,1) 
}
plot(y,type="l",xlab="time")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)

$\phi_1=1.2>1$ の影響は、$\phi_2=-0.8<0$ で相殺され、定常となるようです。
　

標本値と理論値を比べてみましょう。

mv <- c(mean(y),sd(y))
mv #平均と標準偏差の標本値

## [1] 3.228429 2.399325

MV <- c(2/(1-1.2+0.8),sqrt((1+0.8)/(1-0.8)*1^2/((1+0.8)^2-1.2^2)))  
MV #理論値

## [1] 3.333333 2.236068

標本値は理論値と近い値になっています。

3.5.1 自己相関係数・偏自己相関係数

$y_{t}=2+1.2 y_{t-1}-0.8y_{t-2}+\varepsilon _{t}$　のときの理論値は、

$\rho_1=\frac{\phi_1}{1-\phi_2}=\frac{1.2}{1+0.8}=0.667$
$\rho_2=\phi_{1}\frac{\phi_1}{1-\phi_2} + \phi_{2}=\frac{\phi_2+\phi_1^2-\phi_2^2}{1-\phi_2}=\frac{-0.8+1.2^2-(-0.8)^2}{1+0.8}=0$

acf(y,plot=F,lag.max=7)

## 
## Autocorrelations of series 'y', by lag
## 
##      0      1      2      3      4      5      6      7 
##  1.000  0.678  0.018 -0.538 -0.699 -0.447  0.008  0.394

par(family = "HiraKakuProN-W3") #日本語フォントの指定
par(mfrow=c(1,2))
acf(y,plot=T,main="自己相関")
acf(y,plot=T,type="p",main="偏自己相関")

自己相関係数の標本値は理論値とほぼ一致しています。
ラグ１では偏自己相関係数と自己相関係数は同じになります。
間接的影響を与える項が間にないためです。

コレログラムでは、ラグ20までしか表示がありませんが、ラグ200まで表示してみましょう。

par(family = "HiraKakuProN-W3") #日本語フォントの指定
acf(y,plot=T,main="自己相関係数",lag.max=200)

時系列が定常であれば、自己相関係数は 0 に収束していることが分かります。
（逆は成り立ちません）

時系列が定常であれば、誤差項を除いた漸化式は収束するはずです。

　　　$y_{t}=2+1.2 y_{t-1}-0.8y_{t-2}+\varepsilon _{t}$

から誤差項を除いた漸化式をシミュレーションします。

3.5.2 誤差項のない時系列

　　　$y_1=0,~y_2=1,~~~y_{t}=2+1.2 y_{t-1}-0.8y_{t-2}$

y <- rep(0, 200)
y[1] <- 0
y[2] <- 1
for (i in 3:length(y)) {
  y[i] <- 2 + 1.2 * y[i - 1] -  0.8 * y[i - 2]
}

plot(y,type="l",xlab="time")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)

mv <- c(mean(y),sd(y))
mv #平均と標準偏差

## [1] 3.3194444 0.4033081

誤差項のない漸化式で与えられる数列は 0 に収束しています。
変動の様子はコレログラムと似ています。
AR(2) でも、誤差項のない漸化式が収束するとき、時系列は定常になります。

3.5.3 AR(2) の定常性.

3.5.4 非定常AR(2).

　　　$y_1=0,~y_2=1,~~~y_{t}=3 y_{t-1}+10 y_{t-2}+\varepsilon _{t}$

パラメーター $(\phi_1,~\phi_2)=(3,~10)$ は三角形領域の外部にあり、 $y_t$ は非定常です。

y <- rep(0, 120)
y[1] <- 0
y[2] <- 1
for (i in 3:length(y)) {
  y[i] <- 3 * y[i - 1] +  10 * y[i - 2] +rnorm(1,0,1)
}
par(mfrow=c(1,2))
plot(y,type="l",xlab="time",main="Time series")
abline(h=mean(y),col="red")
par(family = "HiraKakuProN-W3") #日本語フォントの指定
acf(y,plot=T,main="ACF",lag.max=120)

時系列は発散しているので定常ではありません。
コレログラムは 0 に収束しているように見えます。
しかし、定常でなければ自己相関係数はラグだけでなく時間にも依存します。
コレログラムは初期値によって変わる可能性があり、コレログラムをもとに統計的判断をすることはできません。

3.5.5 定常AR(2).

　　　$y_1=0,~y_2=1,~~~y_{t}=0.4 y_{t-1}+0.3 y_{t-2}+\varepsilon _{t}$

パラメーター $(\phi_1,~\phi_2)=(0.4,~0.3)$ は三角形領域の内部にあり、 $y_t$ は定常です。

y <- rep(0, 120)
y[1] <- 0
y[2] <- 1
for (i in 3:length(y)) {
  y[i] <- 0.4 * y[i - 1] +  0.3 * y[i - 2]+rnorm(1,0,1)
}
par(mfrow=c(1,2))
plot(y,type="l",xlab="time",main="Time series")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)
par(family = "HiraKakuProN-W3") #日本語フォントの指定
acf(y,plot=T,main="ACF",lag.max=120)

時系列は定常であり、コレログラムは減衰しながら０に収束しているようです。

誤差項を含まない漸化式のグラフと並置すると、

par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(y,type="l",xlab="time",main="誤差項あり")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)

y <- rep(0, 120)
y[1] <- 0
y[2] <- 1
for (i in 3:length(y)) {
  y[i] <- 0.4 * y[i - 1] +  0.3 * y[i - 2]
}
plot(y,type="l",xlab="time",main="誤差項なし")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)

誤差項を除いた漸化式は急速に収束しています。
それに対し、誤差項のある時系列データは、誤差項が影響し大きく変動しています。

参考サイト.
　無粋な日々に様「時系列分析：ARモデルの性質と定常性（1/2）」

定常でない例を、もう一つあげておきます。

3.5.6 非定常AR(2).

　　　$y_1=9,~y_2=5,~~~y_{t}=2+y_{t-1}+0.3y_{t-2}+\varepsilon _{t}$

パラメーター $(\phi_1,~\phi_2)=(1,~0.3)$ は三角形領域の外部にあり、 $y_t$ は非定常です。

y <- rep(0, 120)
y[1] <- 9
y[2] <- 5
for (i in 3:length(y)) {
  y[i] <- 2 + 1 * y[i - 1] +  0.3 * y[i - 2]+rnorm(1,0,1) 
}
par(mfrow=c(1,2))
plot(y,type="l",xlab="time",main="Time series")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)
acf(y,plot=T,main="ACF",lag.max=200)

時系列は発散し、コレログラムは 0 に収束していません。