2017年3月 – かものはしの分析ブログ

はじめに

まわりでシステムトレードや株価の機械学習による予測などに関心が高まってきたので、私も少し調べてみようと思いPythonやRで行われた分析・実装の事例を集めてみました。自分の資産を突っ込む気にはなれないですが、事例を知っておくだけ知っておきたいですね。

調査法

Google検索において以下のクエリで上位に表示されたサイトを中心にまとめました。
「python 機械学習株価」
「python 機械学習為替」
「python アルゴリズム取引」
「python machine learning stock price」
（同様にRも調べました。）

アルゴリズムトレードの理論から学ぶ

システムトレードで億万長者になるぞ! coursera で Computational Investing Part I を受けた
こちらはCourseraで開かれていたシステムトレード向けの講座を受けた方のブログです。Active Portfolio Management: A Quantitative Approach for Producing Superior Returns and Selecting Superior Returns and Controlling Riskはもっとも重要なアルゴリズムトレードの本とされているようです。出版年が1999年と結構古いですが。
Active Portfolio Management: A Quantitative Approach for Producing Superior Returns and Selecting Superior Returns and Controlling Risk (McGraw-Hill Library of Investment and Finance)

Rを用いた株価データのモニタリング

Stock Analysis using R
こちらはアルゴリズム取引とかではないですが、Rで株価を簡単に取得できるquantmodライブラリについて紹介されています。

library(quantmod)
getSymbols("MSFT")
chartSeries(MSFT, subset="last 3 months")
addBBands()

library(quantmod)

getSymbols("MSFT")

chartSeries(MSFT, subset="last 3 months")

addBBands()

以上のコードを実行しただけで時系列データがボックスプロットで描写されました。非常に便利そうです。複数のグラフを並べてディスプレーで見てみたいですね。図ではテクニカル分析で使われるボリンジャーバンドを一発で出してくれています。傾向を掴むための分析に役立つのではないでしょうか。

RSIで株価の連動性を見る

過去のデータからビッグデータ分析で株価を予測する
ここでは、主にpandasを用いて株価の分析を行っています。RSI（Relative Strength Index）という「一定の期間変動幅の中でどれ位株価が上昇しているのか、下落しているのかをはかるもの」を計算して日経平均との相関を見ています。データさえあれば、個別株の市場連動性を見るぶんにはpandasで十分に分析できそうです。

アルゴリズムトレードのシステム構成

自動トレードボット
こちらでは、仮想通貨取引のための自動トレードボットの作成のための手順などが書かれています。
工程としては
・仮想通貨価格データ取得
・バックテストの実施
・明日の価格予想
・学習パラメータの最適化
・結果のメール送信
・APIを使った成行トレード
・ジョブのスケジューリング
などとなっているようです。
完全自動化しようと思うと、作るのは大変そうですね。

WEB屋の自分が機械学習株価予想プログラムを開発した結果
こちらはよりシンプルで
「チャートを見る限り、誰がどうみても今日値上がりする銘柄」を検索して、毎朝Slackに通知してくれるシステム
をPythonで作成されています。最終的に資産を溶かす形になっているようです。

アルゴリズムトレード向けのPythonライブラリ

pythonのアルゴリズムトレードライブラリ
こちらでは
zipline（ http://www.zipline.io/ ）
PyAlgoTrade( http://gbeced.github.io/pyalgotrade/ )
pybacktest( https://github.com/ematvey/pybacktest )
backtrader( https://github.com/mementum/backtrader )
というライブラリが紹介されています。

遺伝的アルゴリズムで為替の自動売買

pythonと遺伝的アルゴリズムで作るFX自動売買システムその1
遺伝的アルゴリズムでFX自動売買その2 進化する売買AIの実装
 遺伝的アルゴリズムでFX自動売買その3 OandaAPIで実際に取引
こちらは、外国為替の取引をPythonで自動化させた試みです。最終的に負けてしまってはいますがシステム構成についても詳しく書かれているので、勉強になりました。データはOANDAというサービスが提供しているAPIを用いて取得し、GA（遺伝的アルゴリズム）を使って為替の売買タイミングを決めているようです。実際の売買にもOandaAPIというものを利用して完全に自動化させています。

決定木による株価リターンの予測

機械学習で未来を予測する – scikit-learn の決定木で未来の株価を予測
こちらはPythonでscikit-learnを用いて、決定木による株価の予測をされています。目的変数としてはリターンインデックスをおいています。前処理においてpandasが使われています。

ランダムフォレストによる株価の予測

Pythonで機械学習を使った株価予測のコードを書こう
こちらでは、ランダムフォレストを用いた機械学習で、ETFなどのデータを特徴量にして個別の株価予測を行っています。

アカデミックサイドでの研究事例

Stock Price Forecasting Using Information from Yahoo Finance and Google Trend
UC　Berkeleyの経済学部生の研究です。こちらはR言語で、ヤフーファイナンスとGoogleトレンドの情報を用いて株価を予測している研究です。これまでの時系列手法よりもパフォーマンスが良いそうです。

最後に

事例を集めてみて、どこまでのレベルのものに手を出すべきか悩ましいと思いました。完全自動化は維持するコストがかかりそうですし、ロスが増大しないか心配です。自分としては情報を自動で取得し、リターンが発生する確率の高そうな銘柄をサジェストしてくれるレベルで十分な気がします。

追記

人工知能学会の全国大会2017でAIを用いた株式運用に関する研究がなされていました。深層学習を用いたものが多そうです。
株価変動パターンの類似性を用いた株価予測
 深層学習と高頻度データを用いた株式注文状況の推定
 LSTM を用いた株価変動予測
 深層学習を用いた株価予測の分析

はじめに

かねてからYahoo!ファイナンスなどの情報に、株価について期間を指定したリターンや標準偏差や回帰分析のあてはまりなどが無いので、そのようなデータを簡単に見れるようにしたいという願望がありました。加えて、できるだけ多くの銘柄の情報をみたいという願望もありました。Rを使えばその全てが簡単に実現できるのでチャレンジしてみました。

・データ収集
・データテーブルの整形
・成長率の計算
・標準偏差の計算
・自己回帰の計算
と非常にシンプルな工程ですが、3600銘柄分のこれらの情報を一気に手に入れれるというのは魅力的だなぁと思います。

データの収集

RFinanceYJというパッケージをインストールして、以下のrfinance_module.Rのスクリプトを置いておきます。

library(RFinanceYJ)

#API
quoteStockTsData &lt;- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval &lt;- substr(time.interval,1,1)
  function.stock &lt;- function(quote.table.item){
    if( xmlSize(quote.table.item) &lt; 5) return(NULL) 
    d &lt;- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    o &lt;- as.number(xmlValue(quote.table.item[[2]]))
    h &lt;- as.number(xmlValue(quote.table.item[[3]]))
    l &lt;- as.number(xmlValue(quote.table.item[[4]]))
    c &lt;- as.number(xmlValue(quote.table.item[[5]]))
    v &lt;- ifelse(xmlSize(quote.table.item) &gt;= 6,as.number(xmlValue(quote.table.item[[6]])),0)
    a &lt;- ifelse(xmlSize(quote.table.item) &gt;= 7,as.number(xmlValue(quote.table.item[[7]])),0)
    return(data.frame(date=d,open=o,high=h,low=l,close=c,volume=v, adj_close=a))
  }
  return(quoteTsData(x,function.stock,since,start.num,date.end,time.interval,type="stock"))
}
quoteFundTsData &lt;- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval &lt;- substr(time.interval,1,1)
  function.fund &lt;- function(quote.table.item){
    d &lt;- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    if(time.interval=='monthly'){
      d &lt;- endOfMonth(d)
    }
    c &lt;- as.number(xmlValue(quote.table.item[[2]]))
    v &lt;- as.number(xmlValue(quote.table.item[[3]]))
    return(data.frame(date=d,constant.value=c,NAV=v))
  }
  return(quoteTsData(x,function.fund,since,start.num,date.end,time.interval,type="fund"))
}
quoteFXTsData &lt;- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval &lt;- substr(time.interval,1,1)
  function.fx &lt;- function(quote.table.item){
    d &lt;- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    o &lt;- as.number(xmlValue(quote.table.item[[2]]))
    h &lt;- as.number(xmlValue(quote.table.item[[3]]))
    l &lt;- as.number(xmlValue(quote.table.item[[4]]))
    c &lt;- as.number(xmlValue(quote.table.item[[5]]))
    return(data.frame(date=d,open=o,high=h,low=l,close=c))
  }
  return(quoteTsData(x,function.fx,since,start.num,date.end,time.interval,type="fx"))
}
######  private functions  #####
#get time series data from Yahoo! Finance.
quoteTsData &lt;- function(x,function.financialproduct,since,start.num,date.end,time.interval,type="stock"){
  r &lt;- NULL
  result.num &lt;- 51
  financial.data &lt;- data.frame(NULL)
  start &lt;- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&amp;sy=\\1&amp;sm=\\2&amp;sd=\\3",since))
  end   &lt;- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&amp;ey=\\1&amp;em=\\2&amp;ed=\\3",date.end))
  
  if(!any(time.interval==c('d','w','m'))) stop("Invalid time.interval value")
  
  extractQuoteTable &lt;- function(r,type){
    if(type %in% c("fund","fx")){
      tbl &lt;- r[[2]][[2]][[7]][[3]][[3]][[9]][[2]]
    }
    else{
      tbl &lt;- r[[2]][[2]][[7]][[3]][[3]][[10]][[2]]
    }
    return(tbl)
  }
  
  #while( result.num &gt;= 51 ){
  while(1){
    start.num &lt;- start.num + 1
    quote.table &lt;- NULL
    quote.url &lt;- paste('http://info.finance.yahoo.co.jp/history/?code=',x,start,end,'&amp;p=',start.num,'&amp;tm=',substr(time.interval,1,1),sep="")
    try(r&lt;-htmlParse(quote.url))
    if( is.null(r) ) stop(paste("Can not access :", quote.url))
    
    try( quote.table &lt;- xpathApply(r,"//table")[[2]], TRUE )
    
    quote.size&lt;-xmlSize(quote.table)
    #cat(paste("size:",quote.size))
    if(xmlSize(quote.table)&lt;=1){
      return (financial.data)
    }
    if( is.null(quote.table) ){
      if( is.null(financial.data) ){
        stop(paste("Can not quote :", x))
      }else{
        financial.data &lt;- financial.data[order(financial.data$date),]
        return(financial.data)
      }
    }
    
    size &lt;- xmlSize(quote.table)
    for(i in 2:size){
      financial.data &lt;- rbind(financial.data,function.financialproduct(quote.table[[i]]))
    }
    
    #result.num &lt;- xmlSize(quote.table)
    Sys.sleep(1)
  }
  financial.data &lt;- financial.data[order(financial.data$date),]
  return(financial.data)  
}
#convert string formart date to POSIXct object
convertToDate &lt;- function(date.string,time.interval)
{
  #data format is different between monthly and dialy or weekly
  if(any(time.interval==c('d','w'))){
    result &lt;- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-\\5",date.string)
  }else if(time.interval=='m'){
    result &lt;- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-01",date.string)
  }
  return(as.POSIXct(result))
}
#convert string to number.
as.number &lt;- function(string)
{
  return(as.double(as.character(gsub("[^0-9.]", "",string))))
}
#return end of month date.
endOfMonth &lt;- function(date.obj)
{
  startOfMonth     &lt;- as.Date(format(date.obj,"%Y%m01"),"%Y%m%d")
  startOfNextMonth &lt;- as.Date(format(startOfMonth+31,"%Y%m01"),"%Y%m%d")
  return(startOfNextMonth-1)
}</pre>

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

library(RFinanceYJ)

#API

quoteStockTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.stock <- function(quote.table.item){

if( xmlSize(quote.table.item) < 5) return(NULL)

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

o <- as.number(xmlValue(quote.table.item[[2]]))

h <- as.number(xmlValue(quote.table.item[[3]]))

l <- as.number(xmlValue(quote.table.item[[4]]))

c <- as.number(xmlValue(quote.table.item[[5]]))

v <- ifelse(xmlSize(quote.table.item) >= 6,as.number(xmlValue(quote.table.item[[6]])),0)

a <- ifelse(xmlSize(quote.table.item) >= 7,as.number(xmlValue(quote.table.item[[7]])),0)

return(data.frame(date=d,open=o,high=h,low=l,close=c,volume=v, adj_close=a))

}

return(quoteTsData(x,function.stock,since,start.num,date.end,time.interval,type="stock"))

}

quoteFundTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.fund <- function(quote.table.item){

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

if(time.interval=='monthly'){

d <- endOfMonth(d)

}

c <- as.number(xmlValue(quote.table.item[[2]]))

v <- as.number(xmlValue(quote.table.item[[3]]))

return(data.frame(date=d,constant.value=c,NAV=v))

}

return(quoteTsData(x,function.fund,since,start.num,date.end,time.interval,type="fund"))

}

quoteFXTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.fx <- function(quote.table.item){

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

o <- as.number(xmlValue(quote.table.item[[2]]))

h <- as.number(xmlValue(quote.table.item[[3]]))

l <- as.number(xmlValue(quote.table.item[[4]]))

c <- as.number(xmlValue(quote.table.item[[5]]))

return(data.frame(date=d,open=o,high=h,low=l,close=c))

}

return(quoteTsData(x,function.fx,since,start.num,date.end,time.interval,type="fx"))

}

###### private functions #####

#get time series data from Yahoo! Finance.

quoteTsData <- function(x,function.financialproduct,since,start.num,date.end,time.interval,type="stock"){

r <- NULL

result.num <- 51

financial.data <- data.frame(NULL)

start <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&sy=\\1&sm=\\2&sd=\\3",since))

end <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&ey=\\1&em=\\2&ed=\\3",date.end))

if(!any(time.interval==c('d','w','m'))) stop("Invalid time.interval value")

extractQuoteTable <- function(r,type){

if(type %in% c("fund","fx")){

tbl <- r[[2]][[2]][[7]][[3]][[3]][[9]][[2]]

}

else{

tbl <- r[[2]][[2]][[7]][[3]][[3]][[10]][[2]]

}

return(tbl)

}

#while( result.num >= 51 ){

while(1){

start.num <- start.num + 1

quote.table <- NULL

quote.url <- paste('http://info.finance.yahoo.co.jp/history/?code=',x,start,end,'&p=',start.num,'&tm=',substr(time.interval,1,1),sep="")

try(r<-htmlParse(quote.url))

if( is.null(r) ) stop(paste("Can not access :", quote.url))

try( quote.table <- xpathApply(r,"//table")[[2]], TRUE )

quote.size<-xmlSize(quote.table)

#cat(paste("size:",quote.size))

if(xmlSize(quote.table)<=1){

return (financial.data)

}

if( is.null(quote.table) ){

if( is.null(financial.data) ){

stop(paste("Can not quote :", x))

}else{

financial.data <- financial.data[order(financial.data$date),]

return(financial.data)

}

size <- xmlSize(quote.table)

for(i in 2:size){

financial.data <- rbind(financial.data,function.financialproduct(quote.table[[i]]))

}

#result.num <- xmlSize(quote.table)

Sys.sleep(1)

}

financial.data <- financial.data[order(financial.data$date),]

return(financial.data)

}

#convert string formart date to POSIXct object

convertToDate <- function(date.string,time.interval)

{

#data format is different between monthly and dialy or weekly

if(any(time.interval==c('d','w'))){

result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-\\5",date.string)

}else if(time.interval=='m'){

result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-01",date.string)

}

return(as.POSIXct(result))

}

#convert string to number.

as.number <- function(string)

{

return(as.double(as.character(gsub("[^0-9.]", "",string))))

}

#return end of month date.

endOfMonth <- function(date.obj)

{

startOfMonth <- as.Date(format(date.obj,"%Y%m01"),"%Y%m%d")

startOfNextMonth <- as.Date(format(startOfMonth+31,"%Y%m01"),"%Y%m%d")

return(startOfNextMonth-1)

}</pre>

以下のコードで2016-01-01から現在に至るまでの、指定した株価コードの株価情報を集めてきます。3600銘柄の株価コードを用意する必要があるので、集めてcodelist.csvで保存しておきます。銘柄コードの情報はここでは掲載しません。（適当なサイトにころがっています。）

source("rfinance_module.R")

codelist <- read.csv("codelist.csv")
since    <- "2016-01-01"

stock_price_database <- 0
for(i in 1:nrow(codelist)){
  if(i ==1){
    dataset <- quoteStockTsData(paste(codelist[i,1],".t",sep = ""),since=since)
    code    <- matrix(codelist[i,1], nrow(dataset), 1)
    dataset <- cbind(code, dataset) 
    stock_price_database <- dataset
  }
  else{
    dataset <- quoteStockTsData(paste(codelist[i,1],".t",sep = ""),since=since)
    code    <- matrix(codelist[i,1], nrow(dataset), 1)
    dataset <- cbind(code, dataset) 
    stock_price_database <- rbind(stock_price_database,dataset)    
  }
}

source("rfinance_module.R")

codelist <- read.csv("codelist.csv")

since <- "2016-01-01"

stock_price_database <- 0

for(i in 1:nrow(codelist)){

if(i ==1){

dataset <- quoteStockTsData(paste(codelist[i,1],".t",sep = ""),since=since)

code <- matrix(codelist[i,1], nrow(dataset), 1)

dataset <- cbind(code, dataset)

stock_price_database <- dataset

}

else{

dataset <- quoteStockTsData(paste(codelist[i,1],".t",sep = ""),since=since)

code <- matrix(codelist[i,1], nrow(dataset), 1)

dataset <- cbind(code, dataset)

stock_price_database <- rbind(stock_price_database,dataset)

}

以下のコードを実行して、元のデータを調整済み終値について、縦が日付で横が銘柄のテーブルに整形します。

#make stock price data table
library(dplyr)
code_list <- data.frame(code=distinct(stock_price_database,code)$code)
date_list <- data.frame(date=distinct(stock_price_database,date)$date)

stock_dataset <- date_list
stock_price_database_adj <- stock_price_database[,which(colnames(stock_price_database) %in% c("code","date","adj_close"))]

for( i in 1:nrow(code_list)){ #nrow(code_list)
  each_stock <- filter(stock_price_database_adj, code == code_list$code[i])
  each_stock <- each_stock[,-1]
  each_stock <- distinct(each_stock, date,.keep_all = TRUE)
  colnames(each_stock)[2] <- paste("sdj",code_list$code[i],sep = "") 
  stock_dataset <- left_join(stock_dataset, each_stock, by = "date")
}

#make stock price data table

library(dplyr)

code_list <- data.frame(code=distinct(stock_price_database,code)$code)

date_list <- data.frame(date=distinct(stock_price_database,date)$date)

stock_dataset <- date_list

stock_price_database_adj <- stock_price_database[,which(colnames(stock_price_database) %in% c("code","date","adj_close"))]

for( i in 1:nrow(code_list)){ #nrow(code_list)

each_stock <- filter(stock_price_database_adj, code == code_list$code[i])

each_stock <- each_stock[,-1]

each_stock <- distinct(each_stock, date,.keep_all = TRUE)

colnames(each_stock)[2] <- paste("sdj",code_list$code[i],sep = "")

stock_dataset <- left_join(stock_dataset, each_stock, by = "date")

}

成長率

指定した期間での成長率を計算します。

#calculate growth rate

#start_date
start_date <- "2016-01-04"
#end_date
end_date <- "2017-02-17"

#return rate datatable
return_data <- filter(stock_dataset, date==start_date |date==end_date)
return_data <- return_data[,-1]
return_data <- rbind(return_data, (as.numeric(return_data[1,]) / as.numeric(return_data[2,])-1)*100)
return_data <- data.frame(code=colnames(return_data),t(return_data))
colnames(return_data) <- c("code",end_date,start_date,"return_rate")
return_data <- return_data[order(return_data$return_rate,decreasing = TRUE),]

#calculate growth rate

#start_date

start_date <- "2016-01-04"

#end_date

end_date <- "2017-02-17"

#return rate datatable

return_data <- filter(stock_dataset, date==start_date |date==end_date)

return_data <- return_data[,-1]

return_data <- rbind(return_data, (as.numeric(return_data[1,]) / as.numeric(return_data[2,])-1)*100)

return_data <- data.frame(code=colnames(return_data),t(return_data))

colnames(return_data) <- c("code",end_date,start_date,"return_rate")

return_data <- return_data[order(return_data$return_rate,decreasing = TRUE),]

標準偏差

ここではシンプルに標準偏差の計算を行います。

#volatility
volatility_data <- stock_dataset[,-1]
volatility_data <- log(volatility_data[-nrow(volatility_data),]) - log(volatility_data[-1,])

volatility <- data.frame(matrix(0,1,2))
colnames(volatility) <- c("code", "stdev")

for(i in 1:length(volatility_data)){
  volatility[i,1] <- colnames(volatility_data)[i]
  volatility[i,2] <- sd(volatility_data[,i],na.rm = TRUE)
}

#volatility

volatility_data <- stock_dataset[,-1]

volatility_data <- log(volatility_data[-nrow(volatility_data),]) - log(volatility_data[-1,])

volatility <- data.frame(matrix(0,1,2))

colnames(volatility) <- c("code", "stdev")

for(i in 1:length(volatility_data)){

volatility[i,1] <- colnames(volatility_data)[i]

volatility[i,2] <- sd(volatility_data[,i],na.rm = TRUE)

}

自己回帰による決定係数

自己回帰でフィッティングの良い銘柄はどれなのかを単純に知りたいので、3600銘柄分計算させます。

#AR(1) estimation

ar1_estimation <- stock_dataset[,-1]
ar1_rsuares <- data.frame(matrix(0,1,2))
colnames(ar1_rsuares) <- c("code", "adj_rsquare")

for(i in 1:length(ar1_estimation)){
  #objective variable
  ar1_y <- ar1_estimation[-nrow(ar1_estimation),i]
  #input variable
  ar1_x <- ar1_estimation[-1,i]
  #combine data
  ar1_data <- data.frame(cbind(ar1_y, ar1_x))
  
  #estimation
  if(class(try(lm(data = ar1_data, formula = ar1_y ~. )))=="lm"){
    model <- lm(data = ar1_data, formula = ar1_y ~. )
    model_sum <- summary(model)
    
    ar1_rsuares[i,1] <- colnames(ar1_estimation)[i]
    ar1_rsuares[i,2] <- model_sum$adj.r.squared    
  }
  else{
    ar1_rsuares[i,1] <- colnames(ar1_estimation)[i]
    ar1_rsuares[i,2] <- "NULL"
  }
}

#AR(1) estimation

ar1_estimation <- stock_dataset[,-1]

ar1_rsuares <- data.frame(matrix(0,1,2))

colnames(ar1_rsuares) <- c("code", "adj_rsquare")

for(i in 1:length(ar1_estimation)){

#objective variable

ar1_y <- ar1_estimation[-nrow(ar1_estimation),i]

#input variable

ar1_x <- ar1_estimation[-1,i]

#combine data

ar1_data <- data.frame(cbind(ar1_y, ar1_x))

#estimation

if(class(try(lm(data = ar1_data, formula = ar1_y ~. )))=="lm"){

model <- lm(data = ar1_data, formula = ar1_y ~. )

model_sum <- summary(model)

ar1_rsuares[i,1] <- colnames(ar1_estimation)[i]

ar1_rsuares[i,2] <- model_sum$adj.r.squared

}

else{

ar1_rsuares[i,1] <- colnames(ar1_estimation)[i]

ar1_rsuares[i,2] <- "NULL"

}

以上より、見てみたかった指標の3600銘柄分を一気に計算することができました。

可視化

試しに一番成長率の高かった株価をプロットしてみます。

#visualization
library(ggplot2)
stock_code <- "sdj7612"

df <- data.frame(stock_dataset[,which(colnames(stock_dataset) %in% c("date", stock_code))])

g <- ggplot(df,aes(x= date, y= df[,2]))+
     geom_line(colour = "magenta",linetype = 2,size = 0.5)+
     geom_smooth (method = "lm")+
     xlab("Date")+
     ylab("Stock Price")+
     ggtitle("Stock Price")
plot(g)

#visualization

library(ggplot2)

stock_code <- "sdj7612"

df <- data.frame(stock_dataset[,which(colnames(stock_dataset) %in% c("date", stock_code))])

g <- ggplot(df,aes(x= date, y= df[,2]))+

geom_line(colour = "magenta",linetype = 2,size = 0.5)+

geom_smooth (method = "lm")+

xlab("Date")+

ylab("Stock Price")+

ggtitle("Stock Price")

plot(g)

おわりに

株価データにカジュアルにアクセスできる環境を手に入れることができたので、今後は様々な指標や時系列解析・業界ごとの分析を試してみたいと思います。

月: 2017年3月

PythonやRを用いたアルゴリズム取引・株価分析まとめ

追記

Rで株価の時系列データを簡単に集計する