Y's note

Web技術・プロダクトマネジメント・そして経営について

本ブログの更新を停止しており、今後は下記Noteに記載していきます。
https://note.com/yutakikuchi/

R言語を用いた自己回帰モデルによる株価予測を試してみた

R 機械学習

一番売れてる株の雑誌ZAiが作った「株」入門改訂版

一番売れてる株の雑誌ZAiが作った「株」入門改訂版

作者: ダイヤモンド・ザイ編集部
出版社/メーカー: ダイヤモンド社
発売日: 2009/03/27
メディア: 単行本
購入: 5人クリック: 71回
この商品を含むブログ (13件) を見る

株価予測

欧州の経済不安により円高/日本株安が深刻になっています。トレーダーとしてはこのBigWaveを見過ごす訳にはいかないですが、「もうはまだなり、まだはもうなり」という言葉があるように投資のタイミングは非常に難しいものです。ここでは投資理論を語るのではなく、機械学習で株価を予測する事を試してみます。今回採用する予測Modelは自己回帰Model(AR)です。ARは時系列データ解析によく用いられます。AR処理はR言語のar関数を用います。

AR(AutoRegressive)Model

ARModel - 自己回帰モデル

ARModelは時系列解析に利用されます。時間と一緒に変動する値に対してARModelを適用し、時系列データに隠れた何かしらの情報を導きだす事を目的としています。AR(p)モデルは次のような式で定義されます。
$X_{t} = c + \sum_{i=1}^{p}\phi_{i} X_{t-i} + \epsilon_{t}$
ここで $X_{t}$ は時刻tで得られた時系列、 $\phi_{i}$ はモデルパラメータ、 $c$ は定数項、 $\epsilon_{t}$ は誤差項となります。定数項は単純化式では削除されることが多いようです。数式から分かるように $X_{t}$ が過去のデータ $X_{t-i}$ に依存していおり、これが回帰と呼ばれる原因になっています。
自己回帰移動平均モデル - Wikipedia

MAModel - 移動平均モデル

MAModelは時系列データの平滑化を行う手法です。誤差による平滑化が行われます。MA(q)は次のような式で定義されます。
$X_{t} = \epsilon_{t} + \sum_{i=1}^{q}\theta_{i}\epsilon_{t-i}$
ここで $X_{t}$ は時刻tで得られた時系列、 $\epsilon_{t}$ は誤差項、 $\theta_{q}$ はモデルパラメータとなります。
移動平均 - Wikipedia

ARMAModel - 自己回帰移動平均モデル

ARMAModelはARModelとMAModelを組み合わせたモデルです。ARMA(p,q)は以下の数式で定義できます。
$X_{t} = \epsilon_{t} + \sum_{i=1}^{q}\theta_{i}\epsilon_{t-i} + \sum_{i=1}^{p}\phi_{i} X_{t-i}$

その他

ここでは詳しく紹介しませんが、その他ARIMA、ARFIMA、ARCH/GARCHなどの時系列解析Modelが存在します。
ARCHモデル - Wikipedia

R言語でのARModel練習

UKgas

R言語に標準で入っているイギリスの1960年〜1986年の四半期毎のガス使用量のデータでARModelの予測練習をしてみます。UKgasのデータは次のものです。また時系列データをplotしてみます。
> UKgas
       Qtr1   Qtr2   Qtr3   Qtr4
1960  160.1  129.7   84.8  120.1
1961  160.1  124.9   84.8  116.9
1962  169.7  140.9   89.7  123.3
1963  187.3  144.1   92.9  120.1
1964  176.1  147.3   89.7  123.3
1965  185.7  155.3   99.3  131.3
1966  200.1  161.7  102.5  136.1
1967  204.9  176.1  112.1  140.9
1968  227.3  195.3  115.3  142.5
1969  244.9  214.5  118.5  153.7
1970  244.9  216.1  188.9  142.5
1971  301.0  196.9  136.1  267.3
1972  317.0  230.5  152.1  336.2
1973  371.4  240.1  158.5  355.4
1974  449.9  286.6  179.3  403.4
1975  491.5  321.8  177.7  409.8
1976  593.9  329.8  176.1  483.5
1977  584.3  395.4  187.3  485.1
1978  669.2  421.0  216.1  509.1
1979  827.7  467.5  209.7  542.7
1980  840.5  414.6  217.7  670.8
1981  848.5  437.0  209.7  701.2
1982  925.3  443.4  214.5  683.6
1983  917.3  515.5  224.1  694.8
1984  989.4  477.1  233.7  730.0
1985 1087.0  534.7  281.8  787.6
1986 1163.9  613.1  347.4  782.8

> ts.plot( UKgas )
ARModel係数

UKgasをARModelに掛けるとパラメータ係数が以下のような結果となります。ar関数はデフォルトではYule-walker方程式により係数を求めるようになっています。
> ar_result <- ar(UKgas)
> ar_result

Call:
ar(x = UKgas)

Coefficients:
      1        2        3        4        5        6  
 0.5284  -0.2951   0.5844   0.3489   0.0436  -0.2709  

Order selected 6  sigma^2 estimated as  8496 
Predict

UKgasの1986年以降のデータを予測させます。predict関数の2番目の引数でどれぐらい先を予測するかを指定できます。ここではn.ahead=40とする事で四半期における10年先の予測をしています。結果の$predが実際の予測値(Prediction)、$seが標準誤差(Standard Error)になります。標準誤差が年数が増えるにつれて少しずつ大きくなっているので、推定精度が悪くなっていることを示しています。
> pred_result <- predict( ar_result, n.ahead = 40 )
> pred_result
$pred
          Qtr1      Qtr2      Qtr3      Qtr4
1987 1053.9970  600.7282  316.9538  748.8121
1988  981.4623  546.8309  308.6325  708.0975
1989  912.3671  504.7132  297.3734  674.2334
1990  849.5541  468.2595  289.8132  643.3941
1991  793.8445  437.2220  284.6219  616.0413
1992  744.0514  411.0465  281.4293  591.6196
1993  699.6330  389.0236  279.8641  569.7661
1994  659.9798  370.6010  279.5886  550.1579
1995  624.5700  355.2798  280.3278  532.5148
1996  592.9406  342.6281  281.8497  516.5963

$se
          Qtr1      Qtr2      Qtr3      Qtr4
1987  92.17482 104.24961 104.25994 111.21667
1988 137.88835 147.29676 147.31935 151.18307
1989 170.10928 175.66459 175.66783 178.40438
1990 191.54920 195.23883 195.24297 197.23555
1991 206.93975 209.38176 209.40298 210.89961
1992 218.25420 219.89611 219.93948 221.09655
1993 226.76889 227.87520 227.93912 228.85094
1994 233.27914 234.02178 234.10233 234.83238
1995 238.31814 238.81219 238.90464 239.49678
1996 242.25722 242.58122 242.68092 243.16651

株価予測

データ取得

Yahoo! Finance - Business Finance, Stock Market, Quotes, News

日経平均株価時系列データ(日足) ichart.finance.yahoo.com

Yahoo! IncのAPIで日経平均のchart情報をcsvで取得することが出来ます。今回は日経平均のデータで予測を試みますが、Yahoo!Incのデータは個別銘柄も取得できるので、個別銘柄でも株価予測は可能だと思います。Yahoo!Incのデータは降順になっているので取得したデータを昇順にsortし直します。
$ wget "http://ichart.finance.yahoo.com/table.csv?s=%5EN225&a=00&b=4&c=1984&d=11&e=31&f=2030&g=d&ignore=.csv" -O kabu.csv
$ head -n 20 kabu.csv
Date,Open,High,Low,Close,Volume,Adj Close
2012-07-27,8543.97,8568.85,8525.19,8566.64,000,8443.10
2012-07-26,8408.22,8448.55,8358.51,8443.10,153400,8443.10
2012-07-25,8410.98,8433.58,8328.02,8365.90,137400,8365.90
2012-07-24,8497.99,8517.52,8446.18,8488.09,126200,8488.09
2012-07-23,8582.28,8606.23,8500.27,8508.32,103400,8508.32
2012-07-20,8785.87,8792.19,8662.72,8669.87,117000,8669.87
2012-07-19,8795.00,8835.80,8771.19,8795.55,113200,8795.55
2012-07-18,8796.15,8802.09,8715.14,8726.74,119000,8726.74
2012-07-17,8740.98,8808.87,8711.73,8755.00,112000,8755.00
2012-07-13,8701.35,8759.06,8695.44,8724.12,122800,8724.12
2012-07-12,8859.39,8862.80,8709.75,8720.01,141800,8720.01
2012-07-11,8819.04,8851.00,8797.73,8851.00,98000,8851.00
2012-07-10,8922.03,8966.99,8855.93,8857.73,108400,8857.73
2012-07-09,8922.99,8954.44,8892.17,8896.88,92000,8896.88
2012-07-06,9052.70,9082.04,8977.35,9020.75,108200,9020.75
2012-07-05,9078.46,9130.54,9069.01,9079.80,97000,9079.80
2012-07-04,9119.69,9136.02,9095.31,9104.17,105600,9104.17
2012-07-03,9013.65,9082.39,9012.86,9066.59,121000,9066.59
2012-07-02,9103.79,9103.79,9003.48,9003.48,99800,9003.48

$ sort kabu.csv | tail -n 143 > kabu_ar.csv
$ cat kabu_ar.csv
Date,Open,High,Low,Close,Volume,Adj Close
2012-01-04,8549.54,8581.45,8547.70,8560.11,106000,8560.11
2012-01-05,8515.66,8519.16,8481.83,8488.71,77600,8488.71
2012-01-06,8488.98,8488.98,8349.33,8390.35,101400,8390.35
2012-01-10,8422.99,8450.59,8405.18,8422.26,112400,8422.26
2012-01-11,8440.96,8463.72,8426.03,8447.88,106200,8447.88
2012-01-12,8423.10,8426.83,8360.33,8385.59,84800,8385.59
2012-01-13,8471.10,8509.76,8458.68,8500.02,109800,8500.02
2012-01-16,8409.79,8409.79,8352.23,8378.36,76400,8378.36
2012-01-17,8420.12,8475.66,8413.22,8466.40,85400,8466.40
2012-01-18,8458.29,8595.78,8446.09,8550.58,136200,8550.58
2012-01-19,8596.68,8668.94,8596.68,8639.68,125000,8639.68
2012-01-20,8751.18,8791.39,8725.32,8766.36,177600,8766.36
2012-01-23,8753.91,8795.27,8744.54,8765.90,119800,8765.90
2012-01-24,8815.36,8825.09,8768.51,8785.33,111000,8785.33
2012-01-25,8842.01,8911.62,8816.09,8883.69,132400,8883.69
2012-01-26,8890.49,8894.60,8834.93,8849.47,122200,8849.47
2012-01-27,8851.02,8886.02,8810.89,8841.22,133200,8841.22
2012-01-30,8803.79,8832.48,8774.23,8793.05,98000,8793.05
予測日数と利用Model

1,3,5,7,10,12,14日後の予測を行い、誤差平均がどれぐらいなのかを見てみます。

今回の予測対象は終値に限定します。

ModelはAR/ARIMAの2つをそれぞれの予測日数に適用します。
学習データと正解データ

学習するデータを2012-01-04〜2012-06-30までを利用します。そこから2012年7月中の株価を予測します。以下は半年間の予測に利用するデータの一部です。これをkabu_ar_training.csvというファイル名で保存します。
(略)
2012-06-01,8465.47,8487.44,8422.50,8440.25,125200,8440.25
2012-06-04,8278.65,8303.35,8238.96,8295.63,126800,8295.63
2012-06-05,8331.02,8388.14,8306.93,8382.00,137000,8382.00
2012-06-06,8428.36,8549.00,8412.55,8533.53,160200,8533.53
2012-06-07,8639.25,8647.79,8599.64,8639.72,127400,8639.72
2012-06-08,8609.78,8611.93,8427.20,8459.26,182800,8459.26
2012-06-11,8612.14,8665.80,8594.58,8624.90,108800,8624.90
2012-06-12,8478.78,8575.87,8452.50,8536.72,113200,8536.72
2012-06-13,8557.57,8615.89,8554.14,8587.84,101600,8587.84
2012-06-14,8531.40,8592.17,8520.99,8568.89,101600,8568.89
2012-06-15,8608.43,8625.19,8552.75,8569.32,106600,8569.32
2012-06-18,8723.55,8766.56,8711.49,8721.02,109800,8721.02
2012-06-19,8692.82,8712.86,8630.66,8655.87,99200,8655.87
2012-06-20,8739.10,8770.41,8711.22,8752.31,114200,8752.31
2012-06-21,8794.03,8859.04,8790.86,8824.07,132600,8824.07
2012-06-22,8733.50,8830.34,8731.79,8798.35,113200,8798.35
2012-06-25,8837.83,8837.83,8726.40,8734.62,92800,8734.62
2012-06-26,8671.61,8712.68,8619.36,8663.99,127600,8663.99
2012-06-27,8677.76,8730.49,8641.52,8730.49,106000,8730.49
2012-06-28,8816.07,8881.44,8805.65,8874.11,118600,8874.11
2012-06-29,8811.46,9044.04,8803.33,9006.78,143800,9006.78
予測しやすいように各データ項目毎(始値、終値、高値、安値)にファイルに落とします。カラムの取得にはunixコマンドのcutを使用します。
$ cut -d ',' -f2 kabu_ar_training.csv > kabu_ar_training_open.txt
$ cut -d ',' -f3 kabu_ar_training.csv > kabu_ar_training_high.txt
$ cut -d ',' -f4 kabu_ar_training.csv > kabu_ar_training_low.txt
$ cut -d ',' -f5 kabu_ar_training.csv > kabu_ar_training_close.txt
予測した結果を正解データと照らし合わせて誤差がどれぐらいかを測定します。以下は2012年7月中の正解データになります。これをkabu_ar_result.csvというファイル名で保存します。正解データも学習データと同じようにカラム毎にファイルを生成しておきます。
$ cat kabu_ar_result.csv
2012-07-02,9103.79,9103.79,9003.48,9003.48,99800,9003.48
2012-07-03,9013.65,9082.39,9012.86,9066.59,121000,9066.59
2012-07-04,9119.69,9136.02,9095.31,9104.17,105600,9104.17
2012-07-05,9078.46,9130.54,9069.01,9079.80,97000,9079.80
2012-07-06,9052.70,9082.04,8977.35,9020.75,108200,9020.75
2012-07-09,8922.99,8954.44,8892.17,8896.88,92000,8896.88
2012-07-10,8922.03,8966.99,8855.93,8857.73,108400,8857.73
2012-07-11,8819.04,8851.00,8797.73,8851.00,98000,8851.00
2012-07-12,8859.39,8862.80,8709.75,8720.01,141800,8720.01
2012-07-13,8701.35,8759.06,8695.44,8724.12,122800,8724.12
2012-07-17,8740.98,8808.87,8711.73,8755.00,112000,8755.00
2012-07-18,8796.15,8802.09,8715.14,8726.74,119000,8726.74
2012-07-19,8795.00,8835.80,8771.19,8795.55,113200,8795.55
2012-07-20,8785.87,8792.19,8662.72,8669.87,117000,8669.87
2012-07-23,8582.28,8606.23,8500.27,8508.32,103400,8508.32
2012-07-24,8497.99,8517.52,8446.18,8488.09,126200,8488.09
2012-07-25,8410.98,8433.58,8328.02,8365.90,137400,8365.90
2012-07-26,8408.22,8448.55,8358.51,8443.10,153400,8443.10
2012-07-27,8543.97,8568.85,8525.19,8566.64,000,8443.10

$ cut -d ',' -f2 kabu_ar_result.csv > kabu_ar_result_open.txt
$ cut -d ',' -f3 kabu_ar_result.csv > kabu_ar_result_high.txt
$ cut -d ',' -f4 kabu_ar_result.csv > kabu_ar_result_low.txt
$ cut -d ',' -f5 kabu_ar_result.csv > kabu_ar_result_close.txt
予測プログラム

AR/ARIMAModelのpredictを書いたRファイルを予め用意します。R上でsource関数にて外部ファイルを読み込みます。Rファイルを以下のように定義します。入力は上で定義した終値の学習データ、出力は予測結果と標準誤差になります。for分で予測したい日付けのpredictを実行しています。少しR言語についても説明をしておくと学習データがtxt形式なのでscan関数で呼び出します。時系列データへの変換はas.ts関数で行います。ar関数でARModel、arima関数でARIMAModelへの適用が可能で、predict関数で何日後まで予測するのかを第二引数で指定します。
kabu_data_close <- as.ts( scan( "/home/yuta/work/R/kabu/data/kabu_ar_training_close.txt" ) )
#ARModel
ar_data_close <- ar( kabu_data_close )
#ARIMAModel
arima_data_close <- arima( kabu_data_close,order=c(1,0,1) )

x = c(1,3,5,7,10,14)
for( i in x ) {
    ar_predict_close <- predict( ar_data_close, n.ahead = i )
    write( ar_predict_close$pred, paste( "/home/yuta/work/R/kabu/data/kabu_ar_predict_close.txt", i ) )
    write( ar_predict_close$se,   paste( "/home/yuta/work/R/kabu/data/kabu_ar_predict_close_se.txt", i ) )

    arima_predict_close <- predict( arima_data_close, n.ahead = i )
    write( arima_predict_close$pred, paste( "/home/yuta/work/R/kabu/data/kabu_arima_predict_close.txt", i ) )
    write( arima_predict_close$se,   paste( "/home/yuta/work/R/kabu/data/kabu_arima_predict_close_se.txt", i ) )
}
データ集計および評価

OpenOfficeに集計したデータをまとめました。ここでは14日後の予測データを上げています。結果から言うと予測値がARもARIMAModelもほぼ線形となっていて正解データと大きくかけ離れてしまっています。時系列は時々刻々と変わるので長期的な予測は厳しいという事が伺えます。時系列の直前のデータを拾わないと意味の無い結果が出てしまうので、再度predictをやり直す事にします。
常に1日後の予測値を求めた評価

上での失敗を活かして常に1日後の予測値を求めます。予測するデータは2012年7月以降、学習データは2012-01-04〜予測の前日までを用います。7月中の14日間繰り返して予測を行い、正解の誤差を比較する事にします。結果を出力するプログラムは次のものを利用します。
kabu_data_close = as.ts( scan( "/home/yuta/work/R/kabu/data/kabu_ar_training_close.txt" ) ) 
ar_predict_result <- matrix()
arima_predict_result <- matrix()
j <- 1
for( i in 123:136) {
    training_kabu_data <- kabu_data_close[1:i]

    #ARModel
    ar_data_close <- ar( training_kabu_data )
    ar_predict_close <- predict( ar_data_close, n.ahead = 1 ) 
    
    #ARIMAModel
    arima_data_close <- arima( training_kabu_data, order=c(1,0,1) )
    arima_predict_close <- predict( arima_data_close, n.ahead = i ) 
    
    ar_predict_result[j] <- ar_predict_close$pred[1]
    arima_predict_result[j] <- arima_predict_close$pred[1]

    j <- j+1 
}
write( ar_predict_result, "/home/yuta/work/R/kabu/data/kabu_ar_predict_close_14.txt" )
write( arima_predict_result, "/home/yuta/work/R/kabu/data/kabu_arima_predict_close_14.txt" )
出力結果をOpenOfficeにまとめました。全体をみると実績値と予測値の波形は同じような形になっています。1日ズレて良い予測値が出ているようです。これは時系列の正解値が与えられるとそこに強く反応して1日後の予測値に影響を与えているようです。AR/ARIMAModel間の差異はそれほど見られませんでした。

まとめ

AR/ARIMAModelを使った株価予測をしてみましたが、様々な予測データを使ってどういった場面に適用可能なのかを追加で調査する必要があると思いました。

14日後などの長期的な予測は線形予測に近づいてしまうため意味の無いものになります。

1日後の予測は当日の正解データに強く影響を受けて、その翌日に近い値が出てしまいます。

1日後の予測でもそれほど当たっているとは言えない結果になっています。

時系列データ解析のModelだけでなく今後NeuralNetworkなどのModelも試して比較評価したいと思いました。

Links

Ｒと時系列(1)

Rの基本パッケージ中の時系列オブジェクト一覧 - RjpWiki

Mase-Rstatman.pdf