www.crosshyou.info

政府統計の総合窓口のデータや、OECDやUCIやのデータを使って、Rの練習をしています。ときどき、読書記録も載せています。

人口の推移と株価の関係の分析1 - R言語でデータの整理(subset関数やmerge関数)

e-Stat(政府統計の総合窓口)のレイアウトが令和になって変更になりました。

f:id:cross_hyou:20190504112145j:plain

今回は右下の「地域」をクリックしてみます。

f:id:cross_hyou:20190504112220j:plain

クリックすると、上図のような画面になります。データ表示をクリックしてみます。

f:id:cross_hyou:20190504112302j:plain

都道府県別の人口を選択して、データを表示します。

f:id:cross_hyou:20190504112343j:plain

このようなデータですね。ダウンロードします。

f:id:cross_hyou:20190504112424j:plain

こんなファイルです。10行目に私が独自の変数名を付けました。

これをread.csvファイルでR言語で読み取りましょう。

f:id:cross_hyou:20190504113230j:plain

株価と人口の関係を調べたいので、株価のファイル、も読込みます。

株価のファイルは下図のようなものです。

f:id:cross_hyou:20190504113522j:plain

これもread.csv関数で読込みます。6行目から読み込むので、skip = 5 としてやります。

f:id:cross_hyou:20190504113840j:plain

株価のデータでは、Timeという変数が西暦を表していて、人口のデータでは、Year1が西暦を表していますが、西暦を百万倍して10万を足した数になっています。これを西暦に直します。

f:id:cross_hyou:20190504114510j:plain

上図のように、うまく変換できました。

次は、merge関数を使って、wdとkabukaを合体させます。

f:id:cross_hyou:20190504115216j:plain

kabukaにYear1という名前で西暦を表す変数を作ってやり、Year1をキーにしてwdとkabukaをmerge関数で統合しました。head関数で始めの数行を表示していますが、1975年がwdは複数行あり、kabukaでは1行だけですが、ちゃんとkabukaのデータがどの1975年にも入っています。

summary関数で基本統計量を表示します。

f:id:cross_hyou:20190504115621j:plain

f:id:cross_hyou:20190504115652j:plain

ここで、分析に必要な変数だけにしておきましょう。必要な変数は、Year1, Pref2, T, M, F, J, JM, JF, KabukaとKabuka1Yです。subset関数を使います。

f:id:cross_hyou:20190504120335j:plain

これでデータの整理は終わりました。

wd, ワーキングデータは、2064の観測と10の変数があります。このうち、KabukaとKabuka1Yがresponse(被説明変数)で、残りがexplanatory(説明変数)です。

1975年から2017年までのデータで、全国と47都道府県のデータです。
Tが総人口、Mが男の総人口、Fが女の総人口、Jが日本人の総人口、JMが男の日本人人口、JFが女の日本人人口です。なので、TとJの差が外国人ということですね。KabukaはTopixの値、Kabuka1Yは1年後の騰落です。1.1だと10%上昇した、ということですね。

これらの人口のデータで株価を説明できるかどうかをこれから調べていきたいと思います。

今回は以上です。

 次回は

 

www.crosshyou.info

 

です。

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

# ファイルの読込み オリジナルデータ = od
od <- read.csv("jinkou.csv", skip = 9, na.strings = c("***", "-", "X"))
str(od)
head(od)
summary(od)

# 株価データの読込み
kabuka <- read.csv("KABUKA_YEARLY.csv", skip = 5)
str(kabuka)
head(kabuka)
summary(kabuka)

# オリジナルデータのYear1を4桁の西暦に直す
wd <- od # オリジナルデータ(od)をワーキングデータ(wd)にコピー
wd$Year1 <- (wd$Year1 - 100000) / 1000000 # 4桁の西暦に変換
head(wd$Year1) # 変換できたか確認

# wdとkabukaの統合
kabuka$Year1 <- kabuka$Time # Timeと同じデータでYear1という変数をkabukaに追加
wd <- merge(wd, kabuka, by = "Year1") # Year1をキーにしてwdとkabukaを統合
str(wd)
head(wd)
summary(wd)

# 必要な変数だけにする
wd <- subset(wd, select = c(Year1, Pref2, T, M, F, J, JM, JF, Kabuka, Kabuka1Y))
str(wd)
head(wd)
summary(wd)

 

読書記録 - 「統計でウソをつく法 - 数式を使わない統計学入門」

 

統計でウソをつく法―数式を使わない統計学入門 (ブルーバックス)

統計でウソをつく法―数式を使わない統計学入門 (ブルーバックス)

 

 私が生まれる前に出版された本です。

平均値が算術平均なのか中央値なのか?という現在ではあまり有効でない記述もありましたが、グラフでウソをつく法などは参考になります、というか昔と全然変わっていないなと思いました。

最終章の統計のウソを見破る五つのカギは心に留めておきたいです。

・誰がそう言っているのか?(統計の出所に注意)

・どういう方法で分かったのか?(調査方法に注意)

・足りないデータはないか?(隠されている資料に注意)

・いっていることが違ってやしないか?(問題のすりかえに注意)

・意味があるかしら?(どこかおかしくないか?)

 

東証一部の規模別・業種別PERとPBRのデータの分析9 - R言語で1社当たりのIncomeとNetAsset

 

www.crosshyou.info

 の続きです。

今回は1社当たりのIncomeとNetAssetの推移を調べてみます。まずは1社当たりのIncome, NetAssetの変数を作成します。

f:id:cross_hyou:20190502162913j:plain

MedianとMeanの値が離れていますね。単位は億円なので、最大の赤字が865億円、最大の黒字が2576億円、平均が216億円、中央値が129億円です。hist関数でヒストグラムを確認します。

f:id:cross_hyou:20190502163132j:plain

f:id:cross_hyou:20190502163143j:plain

だいたい100億円以内のところに集中していますね。

NetAssetはどうでしょうか?

f:id:cross_hyou:20190502163652j:plain

最小は363億円、最大は2兆6264億円、平均値は3244億円、中央値は1988億円です。こちらもhist関数でヒストグラムを描いてみます。

f:id:cross_hyou:20190502164035j:plain

f:id:cross_hyou:20190502164047j:plain

右のすそ野が広い分布ですね。

33業種の中で一番1社当たりの利益が大きい業種はどこでしょうか?

まず、Levelが5だけに絞って、Income1で業種ごとの平均値を出し、小さい順に表示します。

f:id:cross_hyou:20190502165232j:plain

Marine,海運が一番Income1が低く、1社当たり118億円の赤字です。Hoken、保険が一番大きく、1社当たり896億円です。

NetAsset1はどうでしょうか?

f:id:cross_hyou:20190502165617j:plain

Suisan, 農林・水産が一番小さく、1社当たり620億円、Hokenが一番大きく、1社当たり1兆5906億円です。

今回は以上です。

次回は

 

www.crosshyou.info

 

です。

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


# 1社当たりのIncomeの作成
df$Income1 <- df$Income / df$Number
summary(df$Income1)

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

# 1社当たりのNetAssetの作成
df$NetAsset1 <- df$NetAsset / df$Number
summary(df$NetAsset1)

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

# 33業種別のIncme1
wdf <- df[df$Level == 5, ] # 作業用のデータフレーム、wdfを作成
wdf <- tapply(wdf$Income1, wdf$Type, mean, na.rm = TRUE) # 業種別のIncome1の平均値
sort(wdf) # 小さい順に表示
rm(wdf) # wdfを消去

# 33業種別のNetAsset1
wdf <- df[df$Level == 5, ] # 作業用のデータフレーム、wdfを作成
wdf <- tapply(wdf$NetAsset1, wdf$Type, mean, na.rm = TRUE) # 業種別のNetAsset1の平均値
sort(wdf) # 小さい順に表示
rm(wdf) # wdfを消去

東証一部の規模別・業種別PERとPBRのデータの分析8 - R言語で時系列チャートを作成

 

www.crosshyou.info

 の続きです。

今回は、YearとMonthを合わせて、Timeという変数を作ろうと思います。

というのもstr関数でおおもとのデータフレーム、dfを確認すると、

f:id:cross_hyou:20190430164550j:plain

となっていて、YearとMonthが別々になっているので、横軸を時間にして縦軸をPBRとかでグラフを描くのに不便だからです。

Time = Year + (Month - 0.5)/12 という計算式で作成します。

f:id:cross_hyou:20190430165203j:plain

これでX軸をTimeにして、Y軸を様々な変数にしてチャートが作成できます。

まずは、NumberをY軸にしましょう。

f:id:cross_hyou:20190430165952j:plain

f:id:cross_hyou:20190430170008j:plain

col = df$Levelとして、Levelごとに色を分けてグラフにしました。

全般に右肩上がりですね。

次は、PERを描いてみます。

f:id:cross_hyou:20190430170901j:plain

f:id:cross_hyou:20190430170915j:plain

あら。。。Level5だけしか見えないですね。。Level5を除いてグラフにしましょう。

f:id:cross_hyou:20190430171313j:plain

f:id:cross_hyou:20190430171331j:plain

2014年の途中からPERの分布の様子に違いがありますね。

PBRもLevel5を除いてグラフにしてみます。

f:id:cross_hyou:20190430171834j:plain

f:id:cross_hyou:20190430171846j:plain

PBRは14年の途中で分布の様子が変わったということはないですね。

Incomeはどうでしょうか?

f:id:cross_hyou:20190430172209j:plain

f:id:cross_hyou:20190430172226j:plain

Incomeも2014年の途中から水準がバーンと上がっていますね。

NetAssetはどうでしょうか?

f:id:cross_hyou:20190430172524j:plain

f:id:cross_hyou:20190430172538j:plain

こんなグラフになりました。

今回は以上です。
次回は

 

www.crosshyou.info

 

 

です。

 

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


# dfをstr関数で確認
str(df)

# Timeという変数を作成
df$Time <- df$Year + (df$Month - 0.5) / 12
head(df$Time)

# Y軸をNumberにしたチャート
plot(df$Time, df$Number, pch = 21, col = df$Level,
main = "Level 黒:1, 赤:2, 緑:3, 青:4, 水色:5")

# Y軸をPERにしたチャート
plot(df$Time, df$PER, pch = 21, bg = df$Level,
main = "Level 黒:1, 赤:2, 緑:3, 青:4, 水色:5")

# Y軸をPERにしたチャート(Level5を除く)
plot(df$Time[df$Level != 5], df$PER[df$Level != 5],
pch = 21, bg = df$Level[df$Level != 5],
main = "Level 黒:1, 赤:2, 緑:3, 青:4, 水色:5")

# Y軸をPBRにしたチャート(Level5を除く)
plot(df$Time[df$Level != 5], df$PBR[df$Level != 5],
pch = 21, bg = df$Level[df$Level != 5],
main = "Level 黒:1, 赤:2, 緑:3, 青:4, 水色:5")

# Y軸をIncomeにしたチャート(Level5を除く)
plot(df$Time[df$Level != 5], df$Income[df$Level != 5],
pch = 21, bg = df$Level[df$Level != 5],
main = "Level 黒:1, 赤:2, 緑:3, 青:4, 水色:5")

# Y軸をNetAssetにしたチャート(Level5を除く)
plot(df$Time[df$Level != 5], df$NetAsset[df$Level != 5],
pch = 21, bg = df$Level[df$Level != 5],
main = "Level 黒:1, 赤:2, 緑:3, 青:4, 水色:5")

 

東証一部の規模別・業種別PERとPBRデータの分析7 - R言語でヒストグラムを描く

 

www.crosshyou.info

 の続きです。

今回は33業種別のデータで、各変数のヒストグラムを描いてみます。

head関数で基になるデータフレームの始めの15行を表示します。

f:id:cross_hyou:20190430114445j:plain

33業種はLevelが5ですから、まずLevelが5だけのデータフレームを作成します。

f:id:cross_hyou:20190430114734j:plain

まずは、Number、会社数のヒストグラムです。hist関数です。

f:id:cross_hyou:20190430115149j:plain

f:id:cross_hyou:20190430115203j:plain

次は、PERです。

f:id:cross_hyou:20190430115359j:plain

f:id:cross_hyou:20190430115412j:plain

あらら、よくわからないヒストグラムになってしまいましたね。大きな値のPERがあるためでしょう。summary関数で見てみましょう。

f:id:cross_hyou:20190430115620j:plain

最大値で875.70というのがあるからですね。PERが50以下の値だけでヒストグラムを描いてみましょう。

f:id:cross_hyou:20190430115908j:plain

f:id:cross_hyou:20190430115921j:plain

次はPBRのヒストグラムです。

f:id:cross_hyou:20190430120120j:plain

f:id:cross_hyou:20190430120134j:plain

次は、Incomeです。

f:id:cross_hyou:20190430120327j:plain

f:id:cross_hyou:20190430120341j:plain

最後はNetAssetです。

f:id:cross_hyou:20190430120518j:plain

f:id:cross_hyou:20190430120530j:plain

Number, PER, PBR, Income, NetAsset、この5つの変数でバラツキが一番大きい変数は何でしょうか?変動係数を計算してみましょう。変動係数は、標準偏差÷平均値です。

apply関数を使って5つの変数の変動係数をいっぺんに計算しましょう。

f:id:cross_hyou:20190430121856j:plain

こうしてみるとPBRが一番バラツキが小さく、以下、Number、NetAsset、PER、Incomeの順番でバラツキが大きくなっていくことがわかります。

5つの変数間の相関関係はどうでしょうか?cor関数で相関係数を算出してみます。

f:id:cross_hyou:20190430122418j:plain

cor関数を使う前に、na.omit関数でNAのある行を削除しています。

NumberとIncome, NetAssetが以外にも0.5以上の相関係数です。あ、以外でもなんでもないですね。会社数が多ければそれだけ利益や総資産は大きくなりますね。

一社当たりのIncome, NetAssetを計算してみましょう。

f:id:cross_hyou:20190430123033j:plain

この一社当たりのIncome, NetAssetのヒストグラムも描いてみましょう。

まずは、Income1です。

f:id:cross_hyou:20190430123250j:plain

f:id:cross_hyou:20190430123301j:plain

NetAsset1のヒストグラムも描いてみます。

f:id:cross_hyou:20190430123549j:plain

f:id:cross_hyou:20190430123602j:plain

Income1とNetAsset1の変動係数も計算しましょう。

f:id:cross_hyou:20190430125247j:plain

Income1, NetAsset1ともにIncome, NetAssetよりもバラツキは小さくなっています。

変数間の相関係数を見てみましょう。

f:id:cross_hyou:20190430125557j:plain

NumberとIncome1, NetAsset1の相関はマイナス0.2ぐらいとなりました。

今回は以上です。

次回は

 

www.crosshyou.info

 

です。

 

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

# はじめの15行を表示
head(df, 15)

# 33業種だけのデータフレームを作成
df7 <- df[df$Level == 5, ]
head(df7, 15)

# Numberのヒストグラム
hist(df7$Number)

# PERのヒストグラム
hist(df7$PER)

# PERのサマリー
summary(df7$PER)

# 50以下のPERのヒストグラム
hist(df7$PER[df7$PER <= 50])

# PBRのヒストグラム
hist(df7$PBR)

# Incomeのヒストグラム
hist(df7$Income)

# NetAssetのヒストグラム
hist(df7$NetAsset)

# 変動係数の計算
apply(df7[ , c(5:9)], 2, sd, na.rm = TRUE) /
apply(df7[ , c(5:9)], 2, mean, na.rm = TRUE)

# 5つの変数間の相関関係
df8 <- na.omit(df7) #df8は33業種でNAの無いデータフレーム
round(cor(df8[ , c(5:9)]), 3)

# 一社当たりのIncomeとNetAsset
df8$Income1 <- df8$Income / df8$Number
df8$NetAsset1 <- df8$NetAsset / df8$Number

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

# NetAsset1のヒストグラム
hist(df8$NetAsset1)

# Income1の変動係数
sd(df8$Income1) / mean(df8$Income1)

# NetAsset1の変動係数
sd(df8$NetAsset1) / mean(df8$NetAsset1)

# 変数間の相関係数
round(cor(df8[ , c(5:11)]), 3)

 

読書記録 - 「HOPE -- 都市・企業・市民による気候変動総力戦」

 

HOPE――都市・企業・市民による気候変動総力戦

HOPE――都市・企業・市民による気候変動総力戦

  • 作者: マイケル・ブルームバーグ,カール・ポープ
  • 出版社/メーカー: ダイヤモンド社
  • 発売日: 2018/10/18
  • メディア: Kindle版
  • この商品を含むブログを見る
 

 気候変動に対応するには、国家だけでなく都市・企業・市民が自発的に取り組まなくてはならない。そのためには、気候変動を儲かるビジネスにする仕組みを作らなければならない。

この本には世界各国の都市・企業・市民が気候変動問題にいかにして対応しているかの事例が豊富にあるが、残念ながら日本の都市・企業・市民についての言及が一切ない。

それと、マイケルブルームバーグは原子力発電推進派のようだが、福島の原子力発電所の事故とその後の地域住民の生活を知っているのだろうか?

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

 

www.crosshyou.info

 の続きです。

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

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

f:id:cross_hyou:20190427104908j:plain

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

f:id:cross_hyou:20190427105512j:plain

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

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

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

f:id:cross_hyou:20190427110235j:plain

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

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

f:id:cross_hyou:20190427110723j:plain

一番下の行の 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が有意に違うかみてみましょう。

f:id:cross_hyou:20190427111526j:plain

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

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

f:id:cross_hyou:20190427113443j:plain

f:id:cross_hyou:20190427113454j:plain

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

今回は以上です。

 

今回の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")