setwd("D:/R/")

# データファイルを読み込み、
x.tmp <- read.csv("USDJPY-20170831.txt",header=T)
x     <- subset( x.tmp, X.TICKER. == "USDJPY" )$X.CLOSE

# AR(2) モデルを適用する
(fit2 <- arima(x, c(2, 0, 0)))
# 結果中 s.e. は the standard error すなわち標準誤差を表す
# 例えば、http://en.wikipedia.org/wiki/Standard_error_%28statistics%29 が参考になろう
# "intersept" と coefficients ar1 と ar2 の意味は、
# x_n - intersept 
#  = ar1 * (x_{n-1} - intersept) + ar2 * (x_{n-2} - intersept) + epsilon

# (試しに) 時系列（ USDJPY-20170831.txt 中の CLOSE 値）を延長してみよう
# 10 分間（すなわち、1分足を10個）分であるが、AR(2) モデルを用いた予測である。
( predict(fit2, n.ahead=10 ) )

# この結果は「予測」である。
# 真の値を知らないので、この結果を検証する手段がない。
# 従って、この "predict" 関数は、arima関数の評価（実行結果の評価）には使えない
# なぜなら、predict関数は、学習したモデル (上記 fit2 )を入力し, それを学習データ
# の最後の i 個のデータ（ AR(i)モデルの時) に適用し (仮に最後の一個は n番目の
# データであったとしよう）、n+1番目のデータを予測し、この予測結果をあたかも真の
# データとして、最後の i 個のデータにモデルを適用し、n+2番目のデータを予測し、
# それをあたかも真のデータだとして、、、、
# ということを行うため、通常機械学習で行っている予測精度の評価には適用できない。

# というわけで、評価用のプログラムを書くことにした

# USDJPY-20170831.txt をテストデータにする。 

z.tmp <- read.csv("USDJPY-20170831.txt", header=T)
z     <- subset( z.tmp, X.TICKER. == "USDJPY" )$X.CLOSE.

# y.ar を予測値のベクトルとする。予測値は、AR(p) の場合、真のデータの
# p 個を用いて、時間的に次の時点での値を予測する。予測した結果の値は
# 次回以降の予測には用いない。
# 真の値の最初の p 個の次の時点からの予測となる。

# テストデータに対し
y <- z
p <- 5
y.ar <- array( 0, dim = c( length(y) ) )
(fitp <- arima(x, c(p, 0, 0)))
int <- fitp$coef["intercept"]
for ( i in (p+1):length(y) ) { 
  y.ar[i] <- int + coef(fitp)[1:p] %*% (y[(i-1):(i-p)] - int )
}

# データを見ると length(y)=1440
# 予測値 y.ar と真の値 y とを定性的に比較する

tmp <- data.frame(x=y, y=y.ar)
tmp[800:850,]

# 見かけ上、悪くない
# グラフに描いてみよう。真の値 y は黒線、予測値 y.ar は赤線である

plot(10:length(y), y[10:length(y)], type="l")
lines(10:length(y), y[10:length(y)], col=2)

#  y と y.ar はほとんど重なっているようだ。つまり、予測はうまくいっている
# そんなわけは、ない
# もう少し細かく見てみると分る

dev.new()
plotrange= 800:850
plot(plotrange,y[plotrange],type="l")
lines(plotrange,y.ar[plotrange],col=2)

# 明らかに、予測値は、真の値を一個ずらした（さらに少し鈍らした）値になっている。

# テストデータではなく、学習データについても（つまり学習誤差を見ても）同じ
# ことがいえる。
# これは、上記の "y <-z" の部分を "y <- x" としてみれば分る
#
# 次数を変えて、確かめよ。例えば p=1,3,4,5,...