日銀短観2019年12月調査のデータ分析２ - The Central Limit Theorem(中心極限定理)を実感する。

www.crosshyou.info

 の続きです。

今回は、central limit theoremを実感してみたいと思います。日本語だと中心極限定理ですね。

Michael J. Crawley の Statistics: An Introduction Using R

Statistics: An Introduction Using R

には、If you take repeated samples from a population and calculate their averages, then these averages will be normally distributed. とあります。

Now(現状)の分布を見てみます。hist関数でヒストグラムを見てみます。

こんな形状ですね。

こういう分布のデータからデータを1つ抜き出します。これを10000回繰り返します。

こんな感じですね。

それでは、Nowの中から3つサンプルを取り出して、その平均値を計算するのを同じく1万回繰り返してみます。

このように、正規分布のような分布になります。

outcomeは1つの値を取り出したもの、heikinは3つの値を取り出して平均値を計算したものです。

qqnorm関数、qqline関数を実行してみます。

heikinのほうは正規分布とみていい感じですね。

今回は以上です。

日銀短観2019年12月調査のデータ分析１ - R言語でヒストグラムや箱ひげ図を作成する。

今回は先日発表発表された日銀短観2019年12月調査のデータを分析してみようと思います。

短観 : 日本銀行 Bank of Japan

上の画像のZIPファイルの中にあるExcelファイルをダウンロードしました。

こういうファイルです。

これをR言語で分析しようと思います。

こんな感じでCSVファイルにしました。これをread.csv関数でR言語に読込みます。

head関数ではじめの6行だけ表示しました。

summary関数で基本的な統計値を確認しましょう。

Sector, Industry, Sizeがカテゴリカル変数で、Now, NowChg, Next, NextChgが数値変数ですね。短観は0より大きければ景気が良い、0よりも小さければ景気が悪いことを意味します。Nowは現状、Nextは先行きです。どちらの平均値も中央値も0よりも大きいですので、景気は良いのですね。NowChgが現状の変化幅(前回調査との比較)、NextChgが先行きの変化幅です。どちらの平均値、中央値もマイナスですから、前回の調査と比べると景況感は悪化している、ということですね。

数値変数については、分散、標準偏差、変動係数も調べてみましょう。

var関数で分散、sd関数で標準偏差、sd関数 / mean関数で変動係数です。

apply関数と組み合わせて4つの変数についていっぺんに計算しました。

変動係数を見ると、Next(先行き)が一番変動の度合いが大きいことがわかります。

まずは、4つの数値変数について、グラフにしてデータの分布を目でみてみましょう。

まず、上の画像のように、function関数を使って、gという名前の関数を作りました。小さい順のグラフ、箱ひげ図、ヒストグラムの3つのグラフをいっぺんに作成する関数です。

それでは、Now(現状)からみてみましょう。

箱ひげ図を見ると、外れ値は無いことがわかります。

NowChg(現状の変化幅)はどうでしょうか？

上にも下にも外れ値があります。

Next(先行き)はどうでしょうか？

先行きも上下に外れ値がありますね。

NextChg(先行きの変化幅)はどうでしょうか？

上に外れ値がありますね。

４つの変数とも分布の形状は山型で、だいたい左右対称ですね。

今回は以上です。

日経平均とドル円と売買代金の分析７ - 年ごとに日経平均の月次変化に違いはあるのか？棒グラフに標準誤差を表示

www.crosshyou.info

 の続きです。

いままで、YearMonthを無視して分析してきましたので、今回は年によって日経平均の月次変化に違いがあるかどうか、ANOVA分析をしてみたいと思います。

head関数ではじめの数データを表示してみました。左から4文字が年ですね。

この4文字を取り出します。

substr関数ではじめの4文字を取り出します。

2014年はデータが一つですので、分析は無視します。

まず、YearとChgNikkeiだけの作業用のデータフレームを作ります。

NAのある行をna.omit関数で削除したら、2014年はなくなりましたね。tapply関数でYearごとの平均値を出します。

2014年がNAで表示されています。ファクタの水準としてはまだ残っているのですね。整理します。

as.character関数で文字列に戻し、as.factor関数でファクタにまた戻しています。

これで、2014はなくなりました。2018年が0.9939634で一番低く、2019年が1.0097547で一番高いです。

Michael J. CrawleyのStatistics: An Introduction using Rを参考にしてANOVA分析します。

Statistics: An Introduction Using R

 aov関数でANOVA分析ができます。

aov関数でANOVA分析のモデルオブジェクトを作り、summary関数、summary.lm関数で分析結果を表示します。

p値は0.689ですから結論は年ごとの違いがあるとは言えない、となりました。

棒グラフを描いてみます。barplot関数です。

この棒グラフにstandard error(標準誤差)を追加します。

まず、er_barsという関数を作成します。

yが棒グラフのもとのデータ、zが標準誤差、standard errorですね。xは棒グラフのx軸になります。

標準誤差(standard error)はさっきのsummary.lm関数の結果ありましたね。

2016年から2018年は0.015762で2019年は0.016116です。

これで、er_bars関数を使えます、yがheight, zがseです。

上手く標準誤差も棒グラフ上に表示できました。

95%の信頼区間を表示してみましょう。

まず、信頼区間のベクトルを作成します。2015年から2018年はデータ数は12個ですから、自由度は11、2019年はデータ数は11個ですから自由度は10です。

こうして信頼区間のベクトルができたので棒グラフと信頼区間を表示しましょう。

はい。できました。各棒グラフの信頼区間が重なっていますから、各棒グラフに有意な違いは無いことが目で見てわかりますね。

今回は以上です。

日経平均とドル円と売買代金の分析６ - ドル円の水準は日経平均の上げ下げに関連あり。

www.crosshyou.info

 の続きです。

日経平均の前月比変化を前月のドル円や売買代金の水準及び前月比から回帰分析するのは難しいようでした。今回は、前月の日経平均の前月比、前々月のドル円、売買代金の前月比を加えてみます。

まず、NAの行の無い作業用のデータフレームを作成してみました。

lm関数で回帰分析してみます。反応変数はChgNikkeiで説明変数は、PNikkei, PChgNikkei, PYen, PChgYen, PPChgYen, PDaikin, PChgDaikin, PPChgDaikinです。

p-vlueが0.1764と0.05よりも大きいですからダメですね。とりあえず、説明変数を削れるだけ削りましょう。まずは、PChgYenを削ります。

Pr(>F)の列の値が0.7432と0.05よりも大きいですので、model2とmodel1で大きな違いはありません。model2を見てみます。

PPChgYenのp値が一番大きいので、これを削除します。

model3を見てみます。

PDaikinを削ります。

model4を見てみましょう

お！model4のp-valueが0.04855と0.05以下になりました！やっと有意なモデルになりましたね！PChgDaikinを削除します。

model5を見てみます。

PPChgDaikinを削除します。結局前々月の変化幅はモデルには残らないですね。

model6を見てみます。

PNikkeiを削りましょう

model7を見てみます。

PChgNikkeiも削れそうですね。

model8を見てみましょう。

これが一番シンプルなモデルですね。日経平均の前月比は前月のドル円の水準で回帰されます。PYenの係数の符号がマイナスなので、円相場が円安の水準のときは日経平均は下がる、ということですね。

散布図と回帰直線を描きます。

最後にmodel8をグラフ化します。

今回は以上です。

日経平均とドル円と売買代金の分析５ - こんどはロジスティック回帰分析で。こちらもダメ

www.crosshyou.info

 の続きです。

今度はChgNikkeiが1以上なら1、そうでないなら0という2値をとる変数にしてロジスティック回帰分析をしてみます。

まずは2値を取る変数を作成します。

dfChgNikke >= 1 でTRUEとFALSEの論理ベクトルを作り、それをidxとして、作成しました。nichというのが今回の反応変数です。では、やってみましょう。

I(PChgDaikin^2)はいらない感じなので削除します。

Pr(>Chi)の値が0.5934なのでmodel2とmodel1で有意な違いはありません。説明変数の少ないmodel2を採用します。summary関数でみてみましょう。

I(PChgYen^2)は不要のようです。削除したmodel3を作ります。

p値が0.05506なので0.05よりも大きいです。model3を見てみましょう。

う～ん。一番右の列、Pr(>|z|)がみんな0.1以上ですね。。。ロジスティック回帰分析も難しいですね。。。PYen, PDaikin, PNikkeiを入れたものでやってみます。

PChgDaikinはいらないようです。

model2をみてみましょう。

PDaikinはいらない感じですね。

model3をみてみましょう。

PChgYenを削除してみます。

あれ？anova関数でmodel3とmodel4を比較できないですね。。。しばらく考えてわかりました。model3ではNAの行があったけど、model4ではないからですね。NAの無いデータフレームでやり直します。

step関数でいらない変数をいっきに削除しましょう。

最終的に残ったのはPNikkeiとPYenだけでした。

PNikkeiもPYenも係数の符号はマイナスなので、値が大きいほど、nichは0になるということですね。株価が高い・円安のときは株価は下落、株価が低い・円高のときは株価は上昇です。

nichとPNikkei, PYenの散布図を描いてみましょう。

回帰モデルの直線はフィットしてないですね。。。

今回は以上です。

日経平均とドル円と売買代金の分析４ - 前月のドル円と売買代金の変化で今月の日経平均の変化を重回帰分析しようとしたがダメだった。

www.crosshyou.info

 の続きです。

今回は、前月のドル円と売買代金の変化から、今月の日経平均の変化を重回帰分析してみたいと思います。

まず、重回帰分析に必要なデータを用意します。

まず、上の画像のように、PChgYen, PChgDaikinという前月のドル円、売買代金の変化を作成しました。この2つの変数が説明変数で、ChgNikkeiが反応変数です。

まずは散布図を描いてみます。

なんか、全然関係なさそうですね。。

lm関数で重回帰分析をします。まずは、2乗項と交差項も入れたfull modelからスタートします。

p-valueが0.1689と0.05よりも大きいですからダメですね。。。

前月の日経平均、PNikkeiも加えてみましょう。

p-value = 0.1827 とさらに大きくなってしまいました。。全然ダメですね。

PYen, PDaikinも加えて、2乗項は無くしてみましょう。

-- 途中省略 --

p-value = 0.3378 とさらにひどくなりました。

ダメでしたね。。

今回は以上です。

日経平均とドル円と売買代金の分析３ - R言語でヒストグラムや箱ひげ図、quantile - quantile plotを描く

www.crosshyou.info

 の続きです。

今回は月次の変化率のヒストグラムや箱ひげ図を描いてみます。

hist関数でヒストグラムが作成できます。

日経平均は左側、ドル円は右側の裾野が広いですね。売買代金は右側に裾野が広いです。

箱ひげ図はboxplot関数です。

日経平均は小さいほうに外れ値があります。ドル円と売買代金は大きいほうに外れ値があります。

sort関数とplot関数を組み合わせて小さい順に並べてみます。3つ別々では芸がないので、一度にやってみます。

売買代金が一番変動が大きく、ドル円が一番変動が小さいことが一目瞭然ですね。

次は、quantile-quaantile plotを作図してこれらの月次変化が正規分布になっているかどうかをみてみましょう。qqnorm関数とqqline関数を使います。

日経平均とドル円のquantile - quantile plotは両端が点線の直線から外れていますから正規分布では無いことがわかりますね。

shapiro.test関数で正規分布かどうか調べることができます。

shapiro.test関数では、p-valueが0.05より小さいと正規分布とは言えないという結論になります。日経平均と売買代金は正規分布と有意な違いは無いです。それに対してドル円はp-valueが0.007063と0.05よりも小さいので正規分布とは有意に違うということです。日経平均のquantile - qunatile plotぐらいでは正規分布と違うとは言えないのですね。

今回は以上です。