国民経済計算四半期GDP速報データの分析６ - R言語で株価の騰落率とGDPデータの重回帰分析をしてみる。

www.crosshyou.info

 の続きです。

前回は、株価の騰落率とNetExpo(純輸出)で線形単回帰分析をしてみましたが、NetExpoだけでは株価の騰落率を説明することは難しいことがわかりました。そこで今回は、他の変数も加えて、重回帰分析をしてみようと思います。

まずは、株価の騰落率とGDPの各データの相関を調べてみようと思います。cor関数を使います。

NetExpoが-0.3055で相関の絶対値が高かったです。次はMinJuu(民間住宅)が-0.2983ですね。なので、MinJuu(民間住宅)を加えましょう。lm関数を使います。

一番下の行のp-value: 0.001899 がこのモデルの有意かどうかを表しますので、0.05以下なのでモデルは有効だと判断できます。

Interceptのp値は、2.65e-14で文句なく有意、NetExpoは0.0401で0.05以下ですからこれも有意、MinJuuは0.0516ですから5％水準では有意ではないですが、10%水準では有意です。

NetExpo, MinJuuを対数をとったモデルを作ってみましょう。

NetExpoはマイナスの値のデータもありましたので、警告メッセージが出てしまいました。一番下の p-value: 0.05752 を見ると、0.05以上なのでこのモデルは有意ではないですね。

NetExpoとMinJuuの交差項を入れたモデルも調べてみましょう。

一番下の行の p-value: 1.21e-05 なのでモデルは有意です。NetExpo, MinJuu, NetExpo:MinJuuの3つの係数のp値も0.05以下で有意です。

始めのモデル、mlmodel1のAdjusted R-squaredは、0.1118で、交差項も加えたモデル、mlmodel3のAdjusted R-squaredは、0.2249ですから説明力が増しています。

anova関数を使ってmlmodel1とmlmodel3が有意に違うかみてみましょう。

Pr(>F)のところがp値なので、0.0003282と0.05よりも小さい値ですから、2つのモデルに違いはあり、mlmodel3のほうがあてはまりが良いモデルであることがわかりました。

モデルによって予測した株価騰落率と実際の株価騰落率を比較してみましょう。それぞれのモデルの中にfitted.valuesとして格納されていますのでplot関数とlines関数でグラフにしてみます。

グラフにすると、青いライン(交差項のあるモデル)のほうが実際の株価の騰落率の近いことがわかります。

今回は以上です。

今回のR言語のコートは以下にあります。

# Kabuka.1YとGDPデータの相関
cor(data2, data2$Kabuka.1Y)

# Kabuka.1YをNetExpoとMinJuuで重回帰分析
mlmodel1 <- lm(Kabuka.1Y ~ NetExpo + MinJuu, data = data2)
summary(mlmodel1)

# 対数をとったモデル
mlmodel2 <- lm(Kabuka.1Y ~ log(NetExpo) + log(MinJuu), data = data2)
summary(mlmodel2)

# 交差項を入れたモデル
mlmodel3 <- lm(Kabuka.1Y ~ NetExpo * MinJuu, data = data2)
summary(mlmodel3)

# mlmodel1とmlmodel3に有意な違いはあるか
anova(mlmodel1, mlmodel3)

# 株価騰落率とモデルで予測した騰落率
plot(data2$Time, data2$Kabuka.1Y, type = "l",
main = "赤：交差項無し, 青：交差項あり")
lines(data2$Time, mlmodel1$fitted.values, col = "red")
lines(data2$Time, mlmodel3$fitted.values, col = "blue")

国民経済計算四半期GDP速報データの分析５ - R言語で株価、株価の騰落率とGDPデータの相関を調べる。

www.crosshyou.info

 の続きです。

今回は、株価、株価の騰落率とGDPデータの相関係数を調べてみようと思います。

相関係数はcor関数でわかります。

soukan <- の行で、株価とGDPの各データの相関係数を計算しています。round関数で小数点以下３桁までにしています。

names(soukan) <- の行で算出した各データに名前を付けています。

rev(sort の行で相関係数の大きい順に表示しています。Kabukaと一番相関の高いのは、MinSetsu, つまり民間企業設備で、0.526ですね。いちおう確認しましょう。

cor.test関数を使いました。一番下の行の0.5263771とあります。これが相関係数です。さきの計算結果と同じですね。p-value = 1.857e-08なので有意です。

１年後の株価騰落率、Kabuka.1Yとの相関も同じようにやってみましょう。

Kabuka.1YにはNAの行があります。NAがあると相関係数がうまく計算できないので、na.omit関数でNAのないデータフレーム、data2をはじめにつくり、あとは株価との相関係数を算出するのと同じ手順で相関係数を算出しています。絶対値ベースでNetExpo(純輸出)が、-0.306と一番相関がありますね。確認します。

cor.testの結果です。一番下の行の相関係数が -0.305555で一致しています。p-value = 0.003058ですから有意な相関ということです。マイナス相関ですから、純輸出がプラスのときは１年後は株価は安くなっている、純輸出がマイナスのときは１年後は株価は高くなっている、ということですね。

グラフで確認します。散布図を書いてみます。

lm関数で線形単回帰モデルを作り、そのモデルをabline関数で散布図の上に重ねてみました。右肩下がりですが、NetExpoだけで１年後の株価の騰落率を説明するのは難しいようですね。

今回は以上です。

次回は

www.crosshyou.info

です。

今回のR言語のコードです。

# KabukaとGDP各データの相関係数
soukan <- round(cor(data$Kabuka, data[ , 1:12]), 3)
names(soukan) <- colnames(data)[1:12]
rev(sort(soukan))

# KabukaとMinsetsuの相関係数
cor.test(data$Kabuka, data$MinSetsu)

# Kabuka.1YとGDP各データの相関係数
data2 <- na.omit(data)
soukan2 <- round(cor(data2$Kabuka.1Y, data2[ , 1:12]), 3)
names(soukan2) <- colnames(data2)[1:12]
rev(sort(soukan2))

# Kabuka.1YとNetExpoの相関関係
cor.test(data2$Kabuka.1Y, data2$NetExpo)

# Kabuka.1YとNetExpoの散布図
plot(data2$NetExpo, data2$Kabuka.1Y)
abline(lm(Kabuka.1Y ~ NetExpo, data = data2), col = 2)

国民経済計算四半期GDP速報データの分析４ - R言語で各変数のヒストグラムとグラフを作成する。

www.crosshyou.info

 の続きです。

前回まではGDPと株価しかみていませんでしたが、データには、民間最終消費支出や、民間企業設備など様々な変数がありますから、それらを見ていきましょう。

まずは、names関数で変数名を確認します。

MinShouhiは民間最終消費支出です。

hist関数とplot関数でグラフを描きましょう。

民間最終消費支出は山型の分布ですね。

次は、MinJuu(民間住宅)です。でもその前に、ヒストグラムと推移のグラフを作る関数を作ってしまいましょう。function関数です。

こうして、ghp関数というヒストグラムと推移のグラフを描く関数を作成しました。ghp はgraph hist plot という意味です。この関数でMinJuu(民間住宅)を実行します。

うまくできました。民間住宅は右下がりですね。

次は、MinSetsu(民間企業設備)です。

民間企業設備は変動が大きいようです。

MinZaiko(民間在庫変動)はどうでしょうか？

双峰型のヒストグラムです。変動はさらに激しいのかトレンドが無い感じです。

SeiShouhi(政府最終消費支出)を見てみましょう。

政府最終消費支出は右肩上がりのトレンドです。

KoteiShihon(公的固定資本形成)はどうでしょうか？

公的固定資本形成は右肩下がりです。

KouZaiko(公的在庫変動)はどんな感じでしょうか？

公的在庫変動は、2005年ぐらいまでは大きく変動していましたが、以降は変動幅は小さくなっています。

NetExpo(財貨・サービス_純輸出)を見てみましょう。

2010年ぐらいから急激に落ち込んでいる時期があります。

次は、GrossExpo(財貨・サービス_輸出)です。

全体としては右肩上がりですが、2008年ぐらいに急激な落ち込みがありますね。リーマンショックのときでしょう。

GrossImpo(財貨・サービス_輸入)が最後です。

輸入も全体としては右肩上がりですが、輸入と同じように大きな落ち込みがあります。

今回は以上です。

 次回は

www.crosshyou.info

です。

今回のR言語のコードです。

# names関数で各変数の名前を確認
names(data)

# MinShouhi(民間最終消費支出)のグラフ
par(mfrow = c(1, 2))
hist(data$MinShouhi)
plot(data$Time, data$MinShouhi, type = "l")

# ヒストグラムと推移のチャートを描く関数の作成
ghp <- function(x) {
par(mfrow = c(1, 2))
hist(x)
plot(data$Time, x, type = "l")
}

# MinJuu(民間住宅)のグラフ
ghp(data$MinJuu)

# MinSetsu(民間企業設備)のグラフ
ghp(data$MinSetsu)

# MinZaiko(民間在庫変動)のグラフ
ghp(data$MinZaiko)

# SeiShouhi(政府最終消費支出)のグラフ
ghp(data$SeiShouhi)

# KoteiShihon(公的固定資本形成)のグラフ
ghp(data$KoteiShihon)

# KouZaiko(公的在庫変動)のグラフ
ghp(data$KouZaiko)

# NetExpo(財貨・サービス_純輸出)のグラフ
ghp(data$NetExpo)

# GrossExpo(財貨・サービス_輸出)のグラフ
ghp(data$GrossExpo)

# GrossImpo(財貨・サービス_輸入)のグラフ
ghp(data$GrossImpo)

国民経済計算四半期GDP速報データの分析３ - R言語でGDPと株価のチャートをつくる(その２)

www.crosshyou.info

 の続きです。

四半期ごとの株価のデータを景気動向指数の先行系列のデータから取得できたので、その株価データとGDPのデータをくっつけてチャートを作成したいと思います。

このような株価のCSVファイルです。Kbuka-1Yは、1年後にどれだけ株価が上昇しているかを表しでいます。「来年の同じ四半期の株価 / 今年の四半期の株価」で計算しています。

まず、この株価のCSVファイルをread.csv関数で読込みます。skip = 4というオプションで5行目からデータを読込みます。

str関数でデータ構造を、head関数で始めの6行を表示しました。

前回作成したデータフレーム、dataを同じようにstr関数とhead関数で確認しましょう。

Timeの列で一緒にすればいいですね。merge関数を使います。

dataにもkabukaにもKabukaという変数があったので統合後はKabuka.xとKabuka.yになっています。Kabuka.yのほうが正しいほうなので、Kabuka.xを削除します。subset関数を使います。

Kabuka.xがなくなっていますね。

Kabuka.xをKabukaにしましょう。

このように、Kabuka.yをKabukaという名前の新しい変数にデータをコピーして、subset関数でKabuka.yを削除しました。もっと簡単な方法がる気がしますが。

これでようやく四半期の株価がGDPデータに加わりました。plot関数でチャートを作ります。

前回のよりも株価が正確にチャートになりました。

株価をY軸、GDPをX軸に散布図を描き、回帰直線を重ねましょう。

今回は以上です。

 次回は

www.crosshyou.info

です。

今回のR言語のコードです。

# 株価(四半期)の株価を読み込む
kabuka <- read.csv("KABUKA_QUARTERY.csv", skip = 4)
str(kabuka)
head(kabuka)

# dataの確認
str(data)
head(data)

# merge関数でdataとkabukaを統合
data <- merge(data, kabuka, by = "Time")
str(data)
head(data)

# Kabuka.xの削除
data <- subset(data, select = -Kabuka.x)
str(data)

# Kabuka.yをKabukaに変更
data$Kabuka <- data$Kabuka.y
data <- subset(data, select = -Kabuka.y)
str(data)

# GDPと株価のチャート
par(mfrow = c(2,1))
plot(data$Time, data$GDP, col = "red", type = "l")
plot(data$Time, data$Kabuka, col = "blue", type = "l")

# 株価とGDPの散布図と回帰直線
plot(data$GDP, data$Kabuka, pch = 21, col = "red")
abline(lm(data$Kabuka ~ data$GDP), col = "blue")

国民経済計算四半期GDP速報データの分析２ - R言語で株価のデータをくっつけて、GDPと株価のチャートを作る。

www.crosshyou.info

 の続きです。

株価のデータを前回作成したデータフレームにくっつけようと思います。

1994年から2018年末までの株価のデータを探したら、日本取引所グループのウェブサイトに株価の長期時系列のファイルがありました。

こういうファイルです。年初値、年末値、最低値、最高値があるのでこの4つの平均値をその年の平均値にしてしまいましょう。

こんなふうに加工したCSVファイルをR言語のread.csv関数で読み込ませます。

GDPのデータフレームにもこの株価のデータフレームにもTimeという名前の変数があり、これが共通です。このTime変数を基準にして、merge関数で二つのデータフレームを統合します。

できました。

plot関数でGDPとKabukaの相関を見てみましょう。

均等に散らばっている感じですが、cor.test関数で相関関係を調べてみましょう。

相関係数は0.503958でp-valueが9.033e-08ですから正の相関関係で有意ですね。

lm関数で線形単回帰分析をしてみます。

p-valueは9.033e-08なので有意なモデルです。切片、GDPの係数ともにp値は0.05以下なので有意です。R-squareは0.254です。このモデルでは、株価をGDPで25%ぐらい説明できるということですかね。

回帰直線を散布図に重ねてみます。

残差プロットも描きます

横軸をTimeにして推移を描いてみましょう。

GDPはギザギザですね。Kabukaは4つ同じ値が続くので値動きがなめらかというか実際の株価の値動きよりもおだやかです。

今回は以上です。

株価の四半期のデータをどこかからとってきたいですね。

 次回は

www.crosshyou.info

です。

今回のR言語のコードです。

# 株価のCSVファイルの読込み
kabuka <- read.csv("kabuka.csv")
str(kabuka)

# dfとkabukaをTime列で統合
data <- merge(df, kabuka, by = "Time")
head(data)

# GDPとKabukaの散布図
with(data,
plot(GDP, Kabuka)
)

# GDPと株価の相関関係
cor.test(data$GDP, data$Kabuka)

# 株価 = a + b * GDP のモデル
model1 <- lm(Kabuka ~ GDP, data = data)
summary(model1)

# 回帰直線を散布図に重ねる
plot(data$GDP, data$Kabuka, pch = 21, bg = "red")
abline(model1, col = "blue")

# 残差プロット
plot(model1, which = 1)

# 横軸をTimeにしてGDPと株価の推移
par(mfrow = c(2, 1))
plot(data$Time, data$GDP, type = "l", col = "blue")
plot(data$Time, data$Kabuka, type = "l", col = "red")

国民経済計算四半期GDP速報データの分析１ - R言語のread.csv関数でデータを読み込み、分析しやすいようにデータを整える。

今回は、国民経済計算四半期GDP速報のデータをR言語のread.csv関数で読み込んでみたいと思います。

https://www.e-stat.go.jp/dbview?sid=0003109741

このサイトからファイルをダウンロードしました。

こういうファイルです。13行目だけ私が英語の変数名を付け足しています。

このファイルをread.csv関数で読込みます。

read.csv("ファイル名.csv")だと1行目から読み込まれてしまうので、skip = 12というオプションを付けて13行目から読込みます。

また、＊＊＊や - や　ーは数字が得られないものということなので、na.strings = c("＊＊＊", "-", "ー")というオプションを加えてこれらの文字が入っているセルはNAになるように読込みます。

このように始めの12行をスキップして、13行目から読み込むことができました。

Skip1とSkip2の変数はいらないので削除します。

Timeをもう少し簡単にしたいですね。1994000103は1994年の1月から3月期、つまり1994年の第1四半期ですから、1994年のQ1ですね。数字にすると、1994.125が1994年Q1で、1994.375が1994年Q2で、1994.625が1994年Q3で1994.875が1994年Q4です。

とりあえず、このTimeからYearを作成します。

お、うまくいきました。あとは、0.125, 0.375, 0.625, 0.875を足していけばいいですね。R言語はこういうのは簡単です。

R言語は違う長さのベクトルの足し算は、短いベクトルを繰り返して長いベクトルに合わせて計算してくれますからね。

こうして作成したYearをdf$Timeと置き換えてしまいましょう。

TimeJはもういらないので、削除してしまいましょう。

はい、できました。

今回は以上です。

 次回は

www.crosshyou.info

です。

今回のR言語のコードです。

# データの読込み
df <- read.csv("gdp.csv", skip = 12, na.strings = c("***", "-", "－"))
head(df)

# Skip1とSkip2の削除
df <- df[ , -c(1, 2)]
head(df)

# Yearの作成
Year <- round(df$Time / 1000000, 0)
head(Year)

# Yearに小数点以下を足し算
Year <- Year + c(0.125, 0.375, 0.625, 0.875)
head(Year)
tail(Year)

# TimeをYearと置き換え
df$Time <- Year
head(df)

# TimeJを削除
df <- df[ , -2]
head(df)

商業統計調査データの分析５ - R言語のstep関数で回帰モデルの最適な変数を選ぶ。

www.crosshyou.info

 の続きです。

まずは、str関数でデータを確認しましょう。

まず、Yearが何年からあるか確認します。table関数を使います。

1988年から始まって、3年ごとに2012年までの7つの年のデータがあります。

一つの年で、77のデータがあるということは、Sectorが77種類あるということですね。

今回は、Corp_TotalをYear, Staff, Revenue, Inventory, Spaceで回帰分析してみて、どの変数がCorp_Totalに一番影響を与えているのかみてみましょう。

線形回帰分析のlm関数でやってみます。

どの係数も有意でよくわからないですね。

step関数で効果の低い変数を段階的に削除できるらしいです。やってみます。

結局、すべての変数が入っているモデルが一番いいということなのですね。

step関数では、何も変数のないモデルからはじめて、一つずつ変数を追加するという方法もあります。これをやってみます。

まずは変数のないゼロモデルを作成します。

切片は、86112とありますが、これはCorp_Totalの平均値です。

さて、このmodel_zeroからstep関数で変数を追加していきます。

前進ステップワイズ回帰でもやっぱり全部の変数が含まれますね。

今回は以上です。

今回のR言語のコードです。

# データの確認
str(df)

# Yearの確認
table(df$Year)

# Corp_TotalをYear, Staff, Revenue, Inventory, Spaceで回帰
model1 <- lm(Corp_Total ~ Year + Staff + Revenue + Inventory + Space, data = df)
summary(model1)

# step関数
reduced_model <- step(model1, direction = "backward")

# ゼロモデルの作成
model_zero <- lm(Corp_Total ~ 1, data = df)
summary(model_zero)

# Corp_Totalの平均値
mean(df$Corp_Total)

# step関数(前進ステップワイズ回帰)
model_zenshin <- step(model_zero, direction = "forward",
scope = (~ Year + Space + Inventory + Revenue + Staff),
trace = 0)
summary(model_zenshin)