［R］状態空間時系列分析入門 p.2

状態空間時系列分析入門（コマンダー・クープマン著、和合訳）のp.2の計算の再現。

最初に、諸元情報。

> n <- length(UKDriverDeaths)
> x <- 1:n
> y <- log(as.numeric(ts(UKDriverDeaths)))
> k <- 1

切片a、傾き（回帰係数）b、ｙの計算値ye、残差分散（誤差分散）s2、決定係数の二乗r2、F値fval、p値pをそれぞれ求めている。

> g <- matrix(c(rep(1, n), x), n)
> m <- solve(t(g) %*% g) %*% (t(g) %*% y)
> m
             [,1]
[1,]  7.545842732
[2,] -0.001448032
> a <- m[1, 1]
> b <- m[2, 1]
> ye <- a + b * x
> s2 <- sum((y - ye) ^ 2) / (n - 2)
> s2
[1] 0.02299806
> r2 <- 1 - sum((y - ye) ^ 2) / sum((y - mean(y)) ^ 2)
> r2
[1] 0.2205911
> fval <- r2 / (1 - r2) * (n - k - 1) / k
> fval
[1] 53.77447
> p <- pf(fval, k, n - k - 1, lower.tail = FALSE)
> p
[1] 6.313237e-12

Rに標準で搭載されているlm関数を利用する方法もある。

> r <- lm(y ~ 1 + x)
> summary(r)
Call:
lm(formula = y ~ 1 + x)
Residuals:
     Min       1Q   Median       3Q      Max 
-0.33649 -0.10850 -0.01256  0.10368  0.40749 
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.5458427  0.0219747 343.387  < 2e-16 ***
x           -0.0014480  0.0001975  -7.333 6.31e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1517 on 190 degrees of freedom
Multiple R-squared:  0.2206,    Adjusted R-squared:  0.2165 
F-statistic: 53.77 on 1 and 190 DF,  p-value: 6.313e-12
> sum(residuals(r) ^ 2) / (n - 2)
[1] 0.02299806
> r2 <- summary(r)$r.squared
> r2
[1] 0.2205911
> fval <- r2 / (1 - r2) * (n - k - 1) / k
> fval
[1] 53.77447
> p <- pf(fval, k, n - k - 1, lower.tail = FALSE)
> p
[1] 6.313237e-12

計算結果を図にしてみる。

> plot(x,y)
> lines(x, ye)

［R］時系列オブジェクトを使う

Rには標準で時系列解析を行うためのパッケージとその関数群、サンプルデータが含まれている。サンプルデータ（UKDriverDeaths）で試してみる。

> UKDriverDeaths
      Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
1969 1687 1508 1507 1385 1632 1511 1559 1630 1579 1653 2152 2148
1970 1752 1765 1717 1558 1575 1520 1805 1800 1719 2008 2242 2478
1971 2030 1655 1693 1623 1805 1746 1795 1926 1619 1992 2233 2192
1972 2080 1768 1835 1569 1976 1853 1965 1689 1778 1976 2397 2654
1973 2097 1963 1677 1941 2003 1813 2012 1912 2084 2080 2118 2150
1974 1608 1503 1548 1382 1731 1798 1779 1887 2004 2077 2092 2051
1975 1577 1356 1652 1382 1519 1421 1442 1543 1656 1561 1905 2199
1976 1473 1655 1407 1395 1530 1309 1526 1327 1627 1748 1958 2274
1977 1648 1401 1411 1403 1394 1520 1528 1643 1515 1685 2000 2215
1978 1956 1462 1563 1459 1446 1622 1657 1638 1643 1683 2050 2262
1979 1813 1445 1762 1461 1556 1431 1427 1554 1645 1653 2016 2207
1980 1665 1361 1506 1360 1453 1522 1460 1552 1548 1827 1737 1941
1981 1474 1458 1542 1404 1522 1385 1641 1510 1681 1938 1868 1726
1982 1456 1445 1456 1365 1487 1558 1488 1684 1594 1850 1998 2079
1983 1494 1057 1218 1168 1236 1076 1174 1139 1427 1487 1483 1513
1984 1357 1165 1282 1110 1297 1185 1222 1284 1444 1575 1737 1763

tsp関数は、その時系列オブジェクトの開始時刻、終了時刻、観測頻度を返す。

> tsp(UKDriverDeaths)
[1] 1969.000 1984.917   12.000

時刻の基準は時刻の1という量を何に置くかで決まり、1という量を1年でも1時間でも1分としてもどれでもよく、あくまでこれを基準にするということ。このサンプルデータの例では、時刻の基準は年を想定し（1.0=1年）、観測頻度は12としているので、1か月ごとのデータであるということ。つまり、最終データは1984年の12月なので、1984+11/12≒1984.917となる。

start関数とend関数は、それぞれ開始時と終了時の、時系列オブジェクトの基準値と観測頻度における番号を返す。

> start(UKDriverDeaths)
[1] 1969    1
> end(UKDriverDeaths)
[1] 1984   12

開始時データは、基準値が1969で頻度1番目（つまり、1969年1月のデータ）、終了時データは、基準値が1984で頻度12番目（つまり、1984年12月のデータ）ということ。

frequency関数は、その時系列オブジェクトの観測頻度を返す。

> frequency(UKDriverDeaths)
[1] 12

年を基準として毎月のデータであるので、1年は12か月あるので12となっている。

time関数、ts関数は、それぞれ時系列オブジェクトの説明変数、目的変数の値を返す。

> time(UKDriverDeaths)
          Jan      Feb      Mar ・・・
1969 1969.000 1969.083 1969.167 ・・・
1970 1970.000 1970.083 1970.167 ・・・
（表示省略）
> ts(UKDriverDeaths)
Time Series:
Start = 1 
End = 192 
Frequency = 1 
  [1] 1687 1508 1507 1385 1632
  [6] 1511 1559 1630 1579 1653
（表示省略）

as.numeric関数で変換すれば、それぞれベクトルとして値を取り出すことができ、時系列解析ではなく回帰分析も行うことができる。

> x <- as.numeric(time(UKDriverDeaths))
> y <- as.numeric(ts(UKDriverDeaths))
> r <- lm(y ~ x)
> print(r)
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept)            x  
   57107.87       -28.04  
> plot(x, y)
> abline(r)

［R］F分布における分布関数の値を求める

pf関数を使う。以下は、分子の自由度が1、分母の自由度が54の場合の、パーセント点を7.13とした場合の下側（=0.990006…）、上側（=0.009994…）の分布関数の値を求めた例。

> pf(7.13, 1, 54)
[1] 0.9900058
> pf(7.13, 1, 54, lower.tail = FALSE)
[1] 0.009994197

デフォルトでは下側を返す。上側を求めるには、lower.tailオプションをFALSEにする必要がある。

［R］F分布におけるパーセント点を求める

qf関数を使う。以下は、分子の自由度が1、分母の自由度が54の場合の下側1パーセント点（=0.000159…）、上側1パーセント点（=7.128819…）を求めた例。

> qf(0.01, 1, 54)
[1] 0.0001585493
> qf(1 - 0.01, 1, 54)
[1] 7.128819
> qf(0.01, 1, 54, lower.tail = FALSE)
[1] 7.128819

デフォルトでは下側パーセント点を返す。上側パーセント点を求めるには、lower.tailオプションをFALSEとする。