crosshyou

主にクロス表(分割表)分析をしようかなと思いはじめましたが、あまりクロス表の分析はできず。R言語の練習ブログになっています。

全国輸出入コンテナ貨物流動調査のデータの分析 - データファイルを整えるのに苦労する。

いつものように、e-Stat:政府統計の総合窓口を見てみました。

f:id:cross_hyou:20190523201151j:plain

いくつか新着情報がありました。この中で一番マニアックな感じのする、全国輸出入コンテナ貨物流動調査をクリックしてみました。

f:id:cross_hyou:20190523201352j:plain

本調査は、我が国の貿易額の約4割を占める、とあります。コンテナ貨物での貿易って結構需運用難ですね。ファイルをクリックしてみました。

f:id:cross_hyou:20190523201546j:plain

さらにクリックします。

f:id:cross_hyou:20190523201607j:plain

2018年と2013年のデータがあります。2018年を選びました。

f:id:cross_hyou:20190523201747j:plain

たくさん、表がありましたが、この中の表1-4 船積港・船卸港別貨物量 /申告件数/申告価格という表です。

変数名の行を追加しまし不要な行を削除しました。

f:id:cross_hyou:20190523202705j:plain



このファイルをR言語に読込んでデータ分析をしてみます。

まず、read.csv関数で読込みます。

f:id:cross_hyou:20190523202829j:plain

ExPxとInPxが数値型のデータではなく、ファクタ型になってしまっています。

これは問題ですね。

次回に何とかしましょう。

今回は以上です。

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

# データの読込み
df <- read.csv("container.csv")
str(df)

OECDのIndstrial Productionのデータの分析 - R言語での実践例

今回は、OECDのデータを分析してみようと思います。

f:id:cross_hyou:20190522231722j:plain

https://stats.oecd.org/

ここからいろいろデータを取得できるようです。

今回は2017年と2018年のIndustrial Production, 鉱工業生産の伸び率を分析してみます。

f:id:cross_hyou:20190522231945j:plain

Excelにデータをダウンロードしたらこんな感じでした。E8に2017, F8のセルに2018と入力しました。このファイルをR言語のread.csv関数で読み込んで分析しようと思います。

f:id:cross_hyou:20190522232440j:plain

str関数でどんな風に読み込まれたか確認してみます。

f:id:cross_hyou:20190522232552j:plain

うまく読み込んだようです。必要な変数は、Country, X2017, X2018の三つです。subset関数でこの変数だけにしてしまいましょう。

f:id:cross_hyou:20190522232850j:plain

X2017にNAがあるのでNAの行をna.omit関数で削除してしまいます。

f:id:cross_hyou:20190522233035j:plain

X2018が数値データではなくて、ファクタデータになっています。

どんなデータがあるか確認してみます。

f:id:cross_hyou:20190522233255j:plain

空白のものと、.. というものがあるために数値データでなくファクタになってしまっていますね。空白のデータの行を削除します。subset関数です。

f:id:cross_hyou:20190522233710j:plain

空白、.. のデータが0になりました。これを数値データにするには、まず文字列に変換します。as.character関数です。そしてas.numeric関数で数値にします。

f:id:cross_hyou:20190522234041j:plain

X2018が数値データになりました。よく見るとCountryに空白のファクタ水準があります。これを削除します。as.character関数で文字列にしてから、as.factor関数でファクタに戻します。

f:id:cross_hyou:20190522234406j:plain

これでやっとデータが整いました。

2017年のIP(Industrial Production, 鉱工業生産)の平均値は3.332%成長で、2018年は2.105%成長です。この違いは有意な違いなのでしょうか?

まずは、2017年の分散と2018年の分散を比較します。

var.test関数を使います。

f:id:cross_hyou:20190522235011j:plain

p-value = 0.05845 で0.05よりも大きいので分散に違いがあるとは言えないです。

それでは平均値に違いがあるかどうか、t.test関数でみてみましょう。X2017とX2018は同じ国の2017年、2018年のデータがあるので、対になっているデータですね。

f:id:cross_hyou:20190522235519j:plain

p-value = 0.002391 と0.05よりも小さいので、2017年の平均値3.332%と2018年の平均値2.105%は有意な違いがあるということです。

paired = Fとして対のデータとしないで検定してみましょう。

f:id:cross_hyou:20190522235913j:plain

この場合でもp-value = 0.01045と0.05よりも小さいので二つの平均値は明らかに違いがあるとわかります。

視覚的にも違いを見てみましょう。hist関数でヒストグラムを書いてみます。

f:id:cross_hyou:20190523000630j:plain

f:id:cross_hyou:20190523000653j:plain

あきらかに下のヒストグラム、2018年のほうが分布が左によっています。

2017年と2018年の散布図をplot関数で描いてみます。

f:id:cross_hyou:20190523001023j:plain

f:id:cross_hyou:20190523001039j:plain

正の相関があるような気がします。cor.test関数で確認します。

f:id:cross_hyou:20190523001317j:plain


相関係数は0.359でp-valueは0.02111です。p値が0.05より小さいですから、正の相関があるとわかります。

2017年の伸び率の上位の国、下位の国を見てみましょう。order関数を使います。

f:id:cross_hyou:20190523002112j:plain

トルコが9.1%成長で一番です。下位の国はどこでしょうか?

f:id:cross_hyou:20190523002305j:plain

アイルランドが-.2.2%成長で一番低かったです。

2018年の伸び率上位の国はどこでしょうか?

f:id:cross_hyou:20190523002527j:plain

ポーランドが5.9%成長で一番です。あれ?名前が空白の行が残っていましたね。おかしいな。。

最後に2018年の伸び率下位の国はどこでしょうか?

f:id:cross_hyou:20190523002817j:plain

ルクセンブルグが-1.9%成長で一番低いです。

今回は以上です。

R言語のコードは以下のとおりです。

# CSVファイルの読み込み
df <- read.csv("oecd_ip.csv", skip = 7)

# データの確認
str(df)

# 必要な変数だけにする
df <- subset(df, select = c(Country, X2017, X2018))
str(df)
summary(df)

# NA行の削除
df <- na.omit(df)
str(df)
summary(df)

# X2018のデータの確認
table(df$X2018)

# X2018の空白行を削除
df <- subset(df, subset = (X2018 != "") )
df <- subset(df, subset = (X2018 != ".."))
table(df$X2018)

# X2018を文字列に変換
df$X2018 <- as.character(df$X2018)

# X2018を数値に変換
df$X2018 <- as.numeric(df$X2018)

# 確認
str(df)
summary(df)

# Countryの空白のファクタ水準を削除
df$Country <- as.character(df$Country)
df$Country <- as.factor(df$Country)
str(df)
summary(df)

# 分散に違いがあるかどうか
var.test(df$X2017, df$X2018)

# 平均値に違いがあるかどうか
t.test(df$X2017, df$X2018, paired = T)

# 対のデータとみなさいで検証
t.test(df$X2017, df$X2018, paired = F)

# ヒストグラムの作成
par(mfrow = c(2, 1))
hist(df$X2017, breaks = seq(-4,12, 2))
hist(df$X2018, breaks = seq(-4,12, 2))
par(mfrow = c(1, 1))

# 散布図
plot(df$X2017, df$X2018)

# 相関の検定
cor.test(df$X2017, df$X2018)

# 2017年の伸び率上位国
head(df[order(df$X2017, decreasing = TRUE), ])

# 2017年の伸び率下位国
head(df[order(df$X2017), ])

# 2018年伸び率上位国
head(df[order(df$X2018, decreasing = TRUE), ])

# 2018年伸び率下位国
head(df[order(df$X2018), ])

東証一部の規模別・業種別PERとPBRのデータの分析11 - 2013年、2014年と比べるとROEは高くなったがPERは低下している。

 

www.crosshyou.info

 の続きです。

まず、ls()関数で変数として今までどんな変数を作成してきたかを確認します。

f:id:cross_hyou:20190521204200j:plain

dfで始まっているのはデータフレームですね。dfからdf8まで9種類あります。

初心に帰ってdf関数がどんなものか、見てみましょう。head関数で始めの数行を表示します。

f:id:cross_hyou:20190521204624j:plain

Income1が一社当たりの純利益、NetAsset1が一社当たりの純資産です。

Typeの列が規模別・業種別で何を示しているかを表しているのですが、Largeの一社当たりの純利益は945、Smallは16ぐらいと60倍ぐらいの格差があります。これでは東証一部はプレミアム感が無いから再編しよう、という意見が出るのも納得ですね。

試しにIncome1のヒストグラムを描いてみます。hist関数です。

f:id:cross_hyou:20190521205319j:plain

0から250ぐらいでしょうか?そこに分布が集中していますね。

logを取ってからヒストグラムを書いてみましょう。でも、マイナスの値があるからこれを0にする処理をしてからlogを取ってヒストグラムにしてみます。

f:id:cross_hyou:20190521205800j:plain

f:id:cross_hyou:20190521205911j:plain

分布が右に偏っただけで、あんまり分布状況を観察しやすくはなっていないですね。

こんどはTypeをTotalだけに絞り込んで分析しましょう。subset関数を使いました。

f:id:cross_hyou:20190521210311j:plain

PERの時系列の変化を見てみます。plot関数を使います。

f:id:cross_hyou:20190521210645j:plain

f:id:cross_hyou:20190521210659j:plain

plot(df_Total$Time, df_Total$PER)と2回、df_Totalと打ち込むのは面倒だったので、with関数を使いました。

PERは低下傾向ですね。

PBRもグラフにしてみます。

f:id:cross_hyou:20190521211256j:plain

f:id:cross_hyou:20190521211307j:plain

こんどは、with関数のかわりにattach関数を使ってdf_Totalをデフォルトのデータフレームに指定してからplot関数を使いました。detach(df_Total)とするまではdf_Totalがデフォルトのデータフレームになります。PBRはPERほどは低下傾向は示していないようです。

PERとIncome1の散布図を見てみましょう。

f:id:cross_hyou:20190521212009j:plain

f:id:cross_hyou:20190521212019j:plain

同じ色は同じ年を表しています。純利益が大きくなるとPERが低下しています。株価がほぼ変わらずに純利益だけ増えているのでPERが下がっている、ということですね。

PERとNetAsset1はどうでしょうか?

f:id:cross_hyou:20190521212429j:plain

f:id:cross_hyou:20190521212441j:plain

これもIncome1と同じようなグラフですね。

Income1とNetAsset1の散布図を見てみます。

f:id:cross_hyou:20190521212720j:plain

f:id:cross_hyou:20190521212731j:plain

強い正の相関がありますね。Income1をNetAsset1で割ったROEを計算してみます。

f:id:cross_hyou:20190521212951j:plain

ROEの時系列を見てみます。

f:id:cross_hyou:20190521213154j:plain

f:id:cross_hyou:20190521213204j:plain

ROEは上昇トレンドです。

ROEとPERの散布図を見てみます。

f:id:cross_hyou:20190521213608j:plain

f:id:cross_hyou:20190521213637j:plain


ROEが高いほどPERが低いですね。

こうしてみると、PERの低下はファンダメンタル云々というより、2013年、2014年の高いPERのが低下してきて、外国のPERと同じぐらいになった、というように思われます。今回は以上です。

 

 

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

# 変数何があるか確認
ls()

# dfの確認
head(df)

# Income1のヒストグラム
hist(df$Income1)

# Income1の対数変換後のヒストグラム
hist(log(df$Income1 - min(df$Income1)))

# TypeがTotalだけにする
df_Total <- subset(df, subset = (Type == "Total") )
head(df_Total)

# PERの時系列の変化
with(df_Total, plot(Time, PER) )

# PBRの時系列変化
attach(df_Total)
plot(Time, PBR)

# PERとIncome1の散布図
plot(Income1, PER, pch = 21, bg=as.numeric(as.factor(Year)))

# PERとNetAsset1の散布図
plot(NetAsset1, PER, pch = 21, bg = as.numeric(as.factor(Year)))

# Income1とNetAsset1の散布図
plot(Income1, NetAsset1, pch = 21, bg = as.numeric(as.factor(Year)))

# ROEの計算
detach(df_Total)
df_Total$ROE <- df_Total$Income1 / df_Total$NetAsset1

# ROEの時系列
attach(df_Total)
plot(Time, ROE)

# ROEとPERの散布図
plot(ROE, PER, pch = 21, bg = as.numeric(as.factor(Year)))

detach(df_Total)

植生自然度のデータ分析1 - 市街地の割合が一番高いのは大阪府、一番低いのは高知県。

今回はe-Stat(政府統計の総合窓口)から取得した、全国の植生自然度のデータを分析しようと思います。

f:id:cross_hyou:20190518091647j:plain

これがe-Statのサイトの画像ですね。このデータをダウンロードすると、

f:id:cross_hyou:20190518091735j:plain

こうなります。9行目は私が挿入した変数名です。

そもそも、植生自然度ってなんでしょうか?

B 自然環境 | 政府統計の総合窓口

ここをクリックすると定義がありました。

f:id:cross_hyou:20190518092009j:plain

人間の影響がどの程度あるかどうかを示す指標のようです。現在は収集中止だそうです。

早速、データをR言語で読み込んで分析(の真似事)をしてみましょう。

f:id:cross_hyou:20190518093024j:plain

str関数でデータの構造を確認します。

f:id:cross_hyou:20190518093138j:plain

144の観測と13の変数です。Skipする変数があるので、実質は12の変数です。

head関数で始めの数行を表示してみます。

f:id:cross_hyou:20190518093358j:plain

北海道を見るとL9が46.5%となっています。L1が一番人間の手が入っていて、L10が一番自然に近いということが推察されます。

f:id:cross_hyou:20190518094032j:plain

f:id:cross_hyou:20190518094047j:plain

と定義にありました。やっぱりそうですね。

まず、Skip行を削除しましょう。subset関数を使いました。

f:id:cross_hyou:20190518094321j:plain

summary関数でデータの要約統計値を表示します。

f:id:cross_hyou:20190518094558j:plain

Yearを見ると、1988年、1992年、1997年の3つの調査年があることがわかります。もう20年以上も前なんですね。年毎の平均値を見てみましょう。aggregate関数を使いましょう。

f:id:cross_hyou:20190518095122j:plain

L1は1988年度は6.85が1997年度は7.23と増加していますが、その他はそれほど変動が無いように感じます。

次は、各都道府県ごとの平均値を見てみましょう。同じくaggregate関数を使います。

f:id:cross_hyou:20190518095819j:plain

L1、市街地、造成地の割合が一番大きい都道府県はどこでしょうか?rev関数とorder関数をつかってデータフレームを並び替えます。

f:id:cross_hyou:20190518100112j:plain

大阪ですね。東京、神奈川、愛知、埼玉、福岡と続きます。

その反対にL1の割合が一番小さい都道府県はどこでしょうか?

f:id:cross_hyou:20190518100521j:plain

高知が一番、市街地が少ないですね。徳島、北海道、山形、宮崎、島根と続きます。

定義を見ると、L6が植生地、L9が自然林、L10が自然草原でこの三つが人間の手が加わってない場所のようです。この3つを合わせた値が一番大きい都道府県はどこでしょうか?

まず、L6+L9+L10を計算します。

f:id:cross_hyou:20190518101826j:plain

それでは、この新しく作成した変数、L6910の大きい都道府県はどこか調べてみましょう。

f:id:cross_hyou:20190518101708j:plain

北海道が一番です。奈良、宮崎、鹿児島、和歌山、高知と続きます。

反対に小さい都道府県はどこでしょうか?

f:id:cross_hyou:20190518102505j:plain

広島が一番です。大阪、香川、岡山、滋賀、神奈川と続きます。東京が入ってないのと広島が一番というのは意外でした。

今回は以上です。

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


# データファイル(shokusei_shizendo.csv)の読込み
df <- read.csv("shokusei_shizendo.csv", skip = 8, na.strings = c("***", "-", "X") )

# str関数でデータ構造の確認
str(df)

# head関数で始めの数行を表示
head(df)

# Skip列を削除
df <- subset(df, select = -Skip)

# summary関数でデータの要約統計値を表示
summary(df)

# aggregate関数で年度ごとの平均値
aggregate(df[ , -c(1,2)], list(年度 = df$Year), mean)

# aggregate関数で各都道府県ごとの平均値
df_Pref <- aggregate(df[ , -c(1, 2)], list(地域 = df$Pref), mean)
head(df_Pref)

# L1の大きい都道府県はどこか?
head(df_Pref[rev(order(df_Pref$L1)), ])

# L1の小さい都道府県はどこか?
head(df_Pref[order(df_Pref$L1), ])

# L6+L9+L10の作成
df_Pref$L6910 <- df_Pref$L6 + df_Pref$L9 + df_Pref$L10

# L6910の大きい都道府県はどこか?
head(df_Pref[rev(order(df_Pref$L6910)), c("地域", "L6910", "L6", "L9", "L10")])

# L6910の小さい都道府県はどこか?
head(df_Pref[order(df_Pref$L6910), c("地域", "L6910", "L6", "L9", "L10")])

 

人口の推移と株価の関係の分析4 - 埼玉県が人口増加率1位、秋田県が減少率1位。各都道府県の人口増加や男性比率、日本人比率の変化を調べる。

 

www.crosshyou.info

 の続きです。

今回は、人口の伸び率上位・下位、男性比率変化幅の上位・下位、日本人比率変化幅の上位・下位の国を調べようと思います。

まずは、作業用のデータがどんなだったか、確認します。R言語のsummary関数を使います。

f:id:cross_hyou:20190511122710j:plain

Yaer1は1975年から2017年まであります。1975年と2017年を比較して、各都道府県別の変化を見たいので、まず。1975年と2017ねんだけのデータにしましょう。

subset関数をつかいます。

f:id:cross_hyou:20190511123122j:plain

次に、男性比率、日本人比率を表す変数を作成しましょう。

f:id:cross_hyou:20190511123533j:plain

こんどは、分析に必要な変数だけにしましょう。Year1, Pref2, T, MT, JTの5つですね。subset関数を使います。

f:id:cross_hyou:20190511123824j:plain

ここからちょっと頭を使わないとダメですね。まず、Pref2を第1優先、Year1を第2優先でデータを並び替えます。order関数を使います。

f:id:cross_hyou:20190511124516j:plain

head関数で始めの10行を表示しました。1975年の愛知県、2017年の愛知県、1975年の愛媛県、2017年の愛媛県と順番に並んでいます。

このデータフレームに1行ずらした変数を追加します。つまり、1行目が2017年の愛知県のデータになっているものです。

f:id:cross_hyou:20190511125435j:plain

これで、Year1が2017の行はいらなくなったので、削除して、Year1の変数も削除しましょう。

f:id:cross_hyou:20190511125747j:plain

人口の伸び率は、T_1 / Tで計算できますね。男性比率の変化幅は、MT_1 - MTで計算でき、日本人比率の変化幅は、JT_1 - JTで計算できます。

f:id:cross_hyou:20190511130402j:plain

subset関数で必要な変数Pref2, T_CHG, MT_CHG, JT_CHGだけにして、summary関数で見てみましょう。

f:id:cross_hyou:20190511130647j:plain

人口が1.5倍になっている県がありますね。このデータには、全国が入っていますから、全国を削除しましょう。

f:id:cross_hyou:20190511131257j:plain

MT_CHGとJT_CHGを100倍してパーセント表示にしましょう。

f:id:cross_hyou:20190511131516j:plain

hist関数でヒストグラムを作ってみましょう。

f:id:cross_hyou:20190511131712j:plain

f:id:cross_hyou:20190511131724j:plain

1.0以下は人口が減っていることを意味します。19の都道府県で人口が減っています。

男性比率の変化のヒストグラムです。

f:id:cross_hyou:20190511132221j:plain

f:id:cross_hyou:20190511132237j:plain

男性比率が低下している都道府県のほうが多いですね。

次は、日本人比率の変化のヒストグラムです。

f:id:cross_hyou:20190511132456j:plain

f:id:cross_hyou:20190511132508j:plain

日本人比率はほとんどの都道府県で減少していますね。

rev関数とorder関数で人口の伸びの高いところ順に並び替えます。

f:id:cross_hyou:20190511132941j:plain

埼玉が一番の伸びですね!「翔んで埼玉」ですね。一番減少したのは秋田です。

男性比率の変化も高い順に並べてみましょう。

f:id:cross_hyou:20190511133538j:plain

福島が男性比率が一番増加し、北海道が一番減少しています。

次は、日本人比率の変化の大きい順です。

f:id:cross_hyou:20190511133806j:plain

大阪だけが日本人比率が増加しています。東京都が一番日本人比率が減少しています。

pairs関数で3つの変数の散布図を作成してみます。

f:id:cross_hyou:20190511134259j:plain

f:id:cross_hyou:20190511134309j:plain

最後に分析用データフレーム、adを削除して終わります。

f:id:cross_hyou:20190511134703j:plain

今回は以上です。

 

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

# 作業用データの確認
summary(wd)

# 分析用データの作成
# 1975年と2017年だけにする
ad <- subset(wd, subset = (Year1 == 1975 | Year1 == 2017))
summary(ad)

# 男性比率(MT)の作成
ad$MT <- ad$M / ad$T

# 日本人比率(JT)の作成
ad$JT <- ad$J / ad$T

# データの確認
summary(ad[, c("MT", "JT")])

# 分析に必要な変数だけにする
ad <- subset(ad, select = c(Year1, Pref2, T, MT, JT))
summary(ad)

# データの並び替え
ad <- ad[order(ad$Pref2, ad$Year1), ]
head(ad, 10)

# 1行ずらしたデータを追加
ad$T_1 <- c(ad$T[-1], NA)
ad$MT_1 <- c(ad$MT[-1], NA)
ad$JT_1 <- c(ad$JT[-1], NA)
head(ad)
tail(ad)

# 2017の行を削除
ad <- subset(ad, subset = (Year1 != 2017))
head(ad)

# 人口の伸びの計算
ad$T_CHG <- ad$T_1 / ad$T

# 男性比率の変化の計算
ad$MT_CHG <- ad$MT_1 - ad$MT

# 日本人率の変化の計算
ad$JT_CHG <- ad$JT_1 - ad$JT

# 必要な変数だけにする
ad <- subset(ad, select = c(Pref2, T_CHG, MT_CHG, JT_CHG))

# データの確認
summary(ad)

# 全国を削除
ad <- subset(ad, subset = (Pref2 != "全国"))
summary(ad)

# MT_CHGとJT_CHGをパーセント表示に
ad$MT_CHG <- 100 * ad$MT_CHG
ad$JT_CHG <- 100 * ad$JT_CHG
summary(ad)

# T_CHGのヒストグラム
hist(ad$T_CHG)

# MT_CHGのヒストグラム
hist(ad$MT_CHG)

# JT_CHGのヒストグラム
hist(ad$JT_CHG)

# 人口の伸びの高い順
ad[rev(order(ad$T_CHG)), ]

# 男性比率変化の高い順
ad[rev(order(ad$MT_CHG)), ]

# 日本人比率変化の高い順
ad[rev(order(ad$JT_CHG)), ]

# pairs関数で散布図
pairs(ad[ , -1], panel = panel.smooth)

# 分析用データフレーム, adの削除
rm(ad)
ad

人口の推移と株価の関係の分析3- 株価と西暦、総人口、男性比率、日本人比率をR言語のlm関数で重回帰分析する。

 

www.crosshyou.info

 の続きです。

今回はメインタイトルどおり、人口の推移と株価の関係を分析してみようと思います。

まずは、全国の総人口と男性比率と日本人比率、この3つの説明変数と株価の関係を分析してみたいと思います。

まず、分析用のデータフレーム、ad(analysis data.frame)を作成します。

f:id:cross_hyou:20190509193835j:plain

subset関数でまずは全国だけのデータにしました。男性比率(M/T), 日本人比率(J/T)を作成しましょう。

f:id:cross_hyou:20190509194402j:plain

必要な変数は、Year1, T, MT, JTとKabukaです。subset関数で必要な変数だけにしてしまいます。

f:id:cross_hyou:20190509194802j:plain

subset関数で必要な変数だけにして、変数の順番も変更しました。

まずは、pairs関数でそれぞれの変数の散布図を描きます。

f:id:cross_hyou:20190509195143j:plain

f:id:cross_hyou:20190509195155j:plain

グラフマトリックスの一番上の行が株価がY軸で各変数がX軸の散布図です。
cor関数で相関係数マトリックスも作成します。

f:id:cross_hyou:20190509195440j:plain

年と男性比率が一番相関係数の絶対値が大きいですね。男性比率は年々低下しています。

それでは、lm関数で被説明変数を株価、残りを説明変数にして回帰分析をしてみます。

まずは、全ての変数と変数の二乗項、交差項も含んだフルモデルから分析します。

f:id:cross_hyou:20190509200026j:plain

step関数を使って、Kabukaに影響のない変数を削除します。

f:id:cross_hyou:20190509200655j:plain

f:id:cross_hyou:20190509200708j:plain

summary関数で結果を表示します。

f:id:cross_hyou:20190509200828j:plain

あれ、Year1:MTやJTは有意ではないですね。update関数でまず、Year1:JTを削除します。

f:id:cross_hyou:20190509201227j:plain

これで、全ての変数が有意なモデルのなりました。残差をプロットしてみます。

plot関数ですね。

f:id:cross_hyou:20190509201438j:plain

f:id:cross_hyou:20190509201453j:plain

最後に実際の株価とモデルから予測された株価をグラフにしてみましょう。

plot関数で実際の株価のグラフを描いて、points関数でそのグラフにモデルから予想された株価を重ねます。

f:id:cross_hyou:20190509202609j:plain

f:id:cross_hyou:20190509202621j:plain

最後に今回の分析用データ(ad)を消去して終わります。

f:id:cross_hyou:20190509202801j:plain

今回は以上です。

次回は

 

www.crosshyou.info

 

です。

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

# 分析用データ(ad)の作成
ad <- wd # wd(作業用の基データ)をadにコピー
ad <- subset(ad, subset = (Pref2 == "全国")) # 全国だけに絞り込む
head(ad) # 始めの数行の表示

# 男性比率と外国人比率の作成
ad$MT <- ad$M / ad$T # M(男性) / T(総人口)
ad$JT <- ad$J / ad$T # J(日本人) / T(総人口)
summary(ad$MT) # 男性比率(MT)のサマリー
summary(ad$JT) # 日本人比率(JT)のサマリー

# 必要な変数だけにする
ad <- subset(ad, select = c(Kabuka, Year1, T, MT, JT)) # 必要な変数だけ
summary(ad) # サマリー表示

# 変数どうしの散布図
pairs(ad, panel = panel.smooth)

# 相関係数マトリックス
cor(ad)

# フルモデルの分析
full_model <- lm(Kabuka ~ Year1 * T * MT * JT +
I(Year1^2) + I(T^2) + I(MT^2) + I(JT^2), data = ad) # モデルを作成
summary(full_model) # モデルのサマリー

# 不要な変数を削除
reduced_model <- step(full_model, direction = "backward")

# サマリー関数
summary(reduced_model)

# Year1:JTの削除
reduced_model2 <- update(reduced_model, ~ . -Year1:JT)
summary(reduced_model2)

# 残差グラフのプロット
plot(reduced_model2, which = 1)

# 実際の株価とモデルによる株価
plot(ad$Year1, ad$Kabuka, pch = 21, bg = "black", type = "b",
xlab = "西暦", ylab = "株価", main = "黒:実績  赤:モデル値",
ylim = c(0, 3000))
points(reduced_model2$model$Year1, reduced_model2$fitted.values,
pch = 21, bg = "red", type = "b")

# 分析用データ(ad)の消去
rm(ad)
ad

 

東証一部の規模別・業種別PERとPBRのデータの分析10 - PERを利益、資産、時間軸で重回帰分析する

 

www.crosshyou.infoの続きです。

今回はPERを利益、資産、時間軸で重回帰分析したいと思います。

参考書籍は

 

Statistics: An Introduction Using R

Statistics: An Introduction Using R

 

 です。

全体のデータだけで調べようと思います。

まずは、str関数、head関数、summary関数でもともとのデータフレームはどんなデータか確認します。

f:id:cross_hyou:20190506144005j:plain

 3000の観測と12の変数があります。

f:id:cross_hyou:20190506144038j:plain

Levelは1から5まであります。

f:id:cross_hyou:20190506144120j:plain

2013年から2019年ですね。

今回は、Level1, Totalだけで考えたいと思います。

まずは、作業用のデータフレームを作成します。

f:id:cross_hyou:20190506144849j:plain

はじめに、pairs関数でそれぞれの変数どうしの散布図を描きます。

f:id:cross_hyou:20190506145121j:plain

f:id:cross_hyou:20190506145133j:plain

cor関数で相関マトリックスを作成します。

f:id:cross_hyou:20190506145327j:plain

それぞれの変数どうしは相関があることがわかります。

まずは、Time, Income1, NetAsset1とそれぞれの交差項と二乗項を考慮したfull modelから分析しましょう。lm関数を使います。

f:id:cross_hyou:20190506145923j:plain

Time:Income1:NetAsset1は有意ではないようなので、削除したモデルを考えましょう。

update関数を使います。

f:id:cross_hyou:20190506150258j:plain

I(Time^2)は有意ではないようです。これを削除したモデルを考えましょう。

f:id:cross_hyou:20190506150526j:plain

I(NetAsset1^2)も有意ではないようです。削除したモデルを考えましょう。

f:id:cross_hyou:20190506150859j:plain

これですべての係数が有意なモデルとなりました。残差プロットのグラフを描いてみます。。

f:id:cross_hyou:20190506151224j:plain

f:id:cross_hyou:20190506151237j:plain

最後の実績のPERとmodel4で計算されたPERをグラフにして比較します。

f:id:cross_hyou:20190506152116j:plain

f:id:cross_hyou:20190506152127j:plain

今回は以上です。

最後のwdfを消去します。rm関数です。

f:id:cross_hyou:20190506152425j:plain

 次回は

 

www.crosshyou.info

 

です。

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

# データフレームの確認
str(df) # 構造の確認

head(df, 10) # はじめの10行

summary(df) # 各変数のサマリー

# 作業用のデータフレームの作成
wdf <- df[df$Level == 1, c("PER", "Time", "Income1", "NetAsset1")]
str(wdf)
head(wdf, 10)
summary(wdf)

# 変数どうしの散布図
pairs(wdf, panel = panel.smooth)

# 相関マトリックス
cor(wdf)

# full model
model1 <- lm(PER ~ Time * Income1 * NetAsset1 + I(Time^2) + I(Income1^2) +
I(NetAsset1^2), data = wdf)
summary(model1)

# Time:Income1:NetAsset1を削除したモデル
model2 <- update(model1, ~ . -Time:Income1:NetAsset1, data = wdf)
summary(model2)

# I(Time^2)を削除したモデル
model3 <- update(model2, ~ . -I(Time^2), data = wdf)
summary(model3)

# I(NetAsset1^2)を削除したモデル
model4 <- update(model3, ~ . -I(NetAsset1^2), data = wdf)
summary(model4)

# model4の残差プロット
par(mfrow = c(2, 2))
plot(model4)
par(mfrow = c(1, 1))

# X軸Time, Y軸PER
plot(wdf$Time, wdf$PER, pch = 21, bg = "red", main = "赤:実績PER, 青:モデルPER")
points(model4$model2, model4$fitted.values, pch = 21, bg = "blue")

# wdfの削除
rm(wdf)