「回帰分析から学ぶ計量経済学」をPythonで写経 ～ 第２章「結果をどう評価するか」②回帰分析の実践
第２章「結果をどう評価するか」
書籍の著者 山澤成康 先生
この記事は、書籍「回帰分析から学ぶ計量経済学」第２章「結果をどう評価するか」の Python写経活動 を取り扱います。
今回は 回帰分析の推定結果の活用、予測、課題（練習問題） に取り組みます！
線形回帰モデルを使ったデータ分析の基礎体力づくりに相当します。
では書籍を開いて回帰分析の旅に出発です🚀
はじめに
書籍「回帰分析から学ぶ計量経済学」のご紹介
このシリーズは書籍「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」（オーム社、「テキスト」と呼びます）の Python 写経です。
テキストは、2023年11月に発売され、副題「Excelで読み解く経済のしくみ」のとおり、主に Excel を用いて、計量経済学を平易に学べる素晴らしい書籍です。
テキストの「はじめに」に著者の先生が執筆の動機を書かれています。
社会人の統計リテラシーの向上をテーマの１つとした科研費プロジェクトの最終年度で、広く社会人に向けてわかりやすい経済分析の本を書きたかったのです。
私にとって計量経済学は高嶺の花ですが、このテキストでさまざまな回帰分析のアプローチを知ることができました。
また、書籍の Excel 処理を Python に置き換える「寄り道写経」の実践を通じて、回帰分析のお気持ちに少し近づけた感じがいたします。
回帰分析に慣れ親しむのに丁度良いレベル感と内容ですので、これはぜひともブログにしたい！と思って現在に至ります。
計量経済学の色を薄め、データ分析の色を濃いめに書いてまいります！
引用表記
この記事は、出典に記載の書籍に掲載された文章と配布データを引用し、適宜、掲載文章・配布データを改変して書いています。
【出典】
「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」
第１版第１刷、著者 山澤成康、オーム社
記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！
第２章 結果をどう評価するか
この記事は第２章の以下の節を取り扱います。
2.10 消費関数の推定とその判断
2.11 回帰分析による予測
課題1 東京の７月の気温の平均、標準偏差
課題2 出生率低下の要因を探れ
課題3 消費関数と投資関数
記事に用いるデータは、オーム社の書籍紹介サイトからダウンロードできる Excel ファイル内のデータをもとにしてCSVファイルを作成し、data フォルダに格納しています。
第２章で用いるライブラリをインポートします。
### インポート
# 数値計算
import numpy as np
import pandas as pd
# 統計処理
import scipy.stats as stats
import statsmodels.api as sm
import statsmodels.formula.api as smf
# 機械学習(線形回帰・OLS)
from sklearn.linear_model import LinearRegression
# 描画
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'
2.10 消費関数の推定とその判断
テキストのデータ「1994年度から2021年度の国民経済計算データ（内閣府）」を読み込みます。
### データの読み込み
df2 = pd.read_csv('./data/02_02_consumption.csv', index_col=0)
print('df2.shape:', df2.shape)
display(df2.head())【実行結果】
データの２変数の関係について、テキストの公式を引用いたします。
$$
実質民間最終消費 = \alpha + \beta 実質GDP + u
$$
目的変数が実質民間最終消費、説明変数が実質GDPの単回帰になっています！
■ 回帰分析の出力結果 p.78
では回帰分析を実行しましょう！
statsmodels の ols（formula API）を利用します。
大きな値の表示がしやすい summary2() を使います。
### 回帰分析の実行 by statsmodels
result = smf.ols(formula='実質民間最終消費 ~ 実質GDP', data=df2).fit()
display(result.summary2())【実行結果】
続いて回帰の分散分析を実行します。
### 分散分析の実行
display(sm.stats.anova_lm(result).round(4))【実行結果】
■ 用語の違い（Excelと計量経済学） p.79
テキストの「用語の違い」表に沿って、statsmodels の結果をまとめてみます。
$$
\begin{array}{l:l:l:r}
\text{Excel} & 統計用語 & \text{statsmodels} & 推定値 \\
\hline \\
重相関R & 相関係数 & - &- \\
重決定R2 & 決定係数(R^2) &\text{R-squared} & 0.886 \\
補正R2 & 自由度調整済R^2 &\text{Adj. R-squared} & 0.881 \\
観測数 & サンプルサイズ &\text{No. Observations} & 28 \\
観測された分散比 & F値 &\text{F-statistic} & 201.8 \\
有意F & F値のP値 &\text{Prob(F-statistic)} & 0.000 \\
\end{array}
$$
回帰式は以下のようになります。
$$
実質民間最終消費 = 30530.9370 + 0.4971 \times実質GDP + u
$$
係数の検定については以下のようになります。
切片$${\hat{\alpha}}$$の$${P}$$値は$${9.9\%}$$であり、帰無仮説「切片＝ゼロ」を$${5\%}$$水準で棄却できません。
傾き$${\hat{\beta}}$$の$${P}$$値は$${0.0\%}$$であり、帰無仮説「係数＝ゼロ」を$${5\%}$$水準で棄却できます。
2.11 回帰分析による予測
次式のように、推定した係数$${\hat{\alpha}, \hat{\beta}}$$を用いて目的変数の予測値（テキストの理論値）$${\hat{Y}_i}$$を計算できます。
$$
\hat{Y}_i = \hat{\alpha} + \hat{\beta} X_i
$$
■ 実績値と理論値 p.81
「1994年度から2021年度の国民経済計算データ（内閣府）」の説明変数：実質GDPを用いて、実質民間最終消費の予測値を推定し、テキストの「図表 実績値と理論値」のチャートを描画します。
statsmodels では、推定結果 result の属性 fittedvalues で予測値（理論値）を取得できます。
### 観測値と推定値（予測値）の描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 4))
# 実績値の折れ線グラフの描画
ax.plot(df2['実質民間最終消費'], label='実績値')
# 理論値の折れ線グラフの描画 ※result.fittedvaluesで学習に用いたXによる予測値を取得
ax.plot(result.fittedvalues, color='tab:red', ls='--', label='理論値')
# 残差の棒グラフの描画 ※result.residで残差を取得
ax.bar(df2.index, result.resid, color='tab:green', alpha=0.8, label='残差')
# 縦軸0の水平線の描画
ax.axhline(0, color='black', lw=0.8)
# 修飾
ax.set(ylim=(-50000, 350000), xlabel='年度', ylabel='実質民間最終消費 [10億円]',
title='実質民間最終消費の実績値・理論値・残差')
ax.grid(lw=0.5)
ax.legend()
plt.show()【実行結果】
また、説明変数を与えて目的変数の予測をする場合、statsmodels では、推定結果（モデル） result の predict メソッドで計算できます。
次のコードは、実質GDPの値を 590000 から 700000 まで 10000 間隔で設定して、実質民間最終消費の予測値を計算しています。
### statsmodelsの線形回帰での予測
# X:実質GDPの値の設定 ※formulaを用いたのでデータフレームと同じ変数名を付ける
X_new = dict(実質GDP=list(range(590000, 700001, 10000)))
# 予測の実行
y_pred = result.predict(exog=X_new)
# データフレーム化して表示
display(pd.DataFrame(
{'実質GDP': X_new['実質GDP'], '実質民間最終消費（予測値）': y_pred}))【実行結果】
課題1 東京の７月の気温の平均、標準偏差
気象庁のWebサイトから1875年から2023年までの東京都の月次平均気温データを取得して、７月の平均気温の平均と標準偏差、2023年7月の$${Z}$$値とパーセンタイル点（正規分布に従うと仮定）を算出します。
第１章の課題１で氏油した月次平均気温CSVデータを読み込みます。
### データの読み込み
df3 = pd.read_csv('./data/01_ex1.csv', index_col=0)
print('df3.shape:', df3.shape)
display(df3.head())【実行結果】
7月の平均気温の平均値・標準偏差を計算します。
### 7月の平均気温の平均値・標準偏差の計算
# 平均値の計算
mean_july = df3['7月'].mean()
# 不偏分散の標準偏差の計算
std_july = df3['7月'].std(ddof=1)
# 結果の表示
print('平均値 : ', mean_july)
print('標準偏差: ', std_july)【実行結果】
平均値は約 25.1℃、標準偏差は約 1.6 ℃です。
続いて、2023年7月の$${Z}$$値とそのパーセンタイル点を計算します。
### 2023年7月の平均気温のZ値と標準正規分布の下位からの％を算出
# z値の計算
z_value = (df3.loc[2023, '7月'] - mean_july) / std_july
# 下位からの％を算出 ※標準正規分布の累積分布関数cdfを使用
z_percent = stats.norm.cdf(z_value, loc=0, scale=1)
# 結果の表示
print('ｚ値 : ', z_value)
print('下位からの％: ', z_percent)【実行結果】
2023年は 98.9%点です。
過去の7月平均気温と比べてずいぶん暑くなっているようです。
課題2 出生率低下の要因を探れ
第２章 2.8 節で使用した合計特殊出生率CSVデータを読み込み、目的変数：合計特殊出生率、説明変数：女性未婚率・有配偶出生率として回帰分析を行います。
### データの読み込み
df4 = pd.read_csv('./data/02_01_Ftest.csv', index_col=0)
df4.columns = ['合計特殊出生率', '女性未婚率', '有配偶出生率', '定数項回帰用']
print('df4.shape:', df4.shape)
display(df4.head())【実行結果】
回帰分析を実行します。
### 回帰分析
# 回帰分析の実行
result = smf.ols(
formula = '合計特殊出生率 ~ 女性未婚率 + 有配偶出生率',
data=df4
).fit()
# 結果の表示
display(result.summary())【実行結果】
次の３点のみに注目すると回帰モデルはいい感じ、でしょうか。
$${F}$$値：$${P}$$値が 0.0 %なので回帰は有意（ゼロでない係数がある）
全ての係数の$${t}$$値：$${p}$$値が 0.0 %なので係数はゼロではない
自由度調整済み決定係数：0.910 でありモデルの当てはまりは良い
係数を取り出してみましょう。
display(result.params.rename('係数').to_frame())【実行結果】
係数の推定値から次のことが推察できそうです。
女性未婚率が 1pt 増加すると、合計特殊出生率は 0.028pt 減少する
有配偶出生率が 1pt 増加すると、合計特殊出生率は 0.014pt 増加する
課題3 消費関数と投資関数
テキストの国民経済計算関連データを読み込んで、消費関数と投資関数の回帰分析を行い、決定係数と$${t}$$値を調べて、どちらの関数が当てはまりがよいか調べます。
### データの読み込み
df5 = pd.read_csv('./data/02_ex1.csv', index_col=0)
print('df5.shape:', df5.shape)
display(df5.head())【実行結果】
以下の消費関数に関して、回帰分析を実行します。
$$
実質民間最終消費 = \alpha + \beta \times 実質GDP + 誤差
$$
### 消費関数の回帰分析 消費＝α＋β×実質GDP
# 回帰分析
result = smf.ols(
formula='実質民間最終消費 ~ 実質GDP', data=df5).fit()
# 結果表示
display(result.summary())【実行結果】
続いて、以下の投資関数に関して、回帰分析を実行します。
$$
実質民間設備投資 = \alpha + \beta \times 国内銀行貸出約定平均金利 + 誤差
$$
### 投資関数の回帰分析 投資＝α＋β×金利
# 回帰分析
result = smf.ols(
formula='実質民間設備投資 ~ 国内銀行貸出約定平均金利', data=df5).fit()
# 結果表示
display(result.summary())【実行結果】
消費関数と投資関数を比較します。
モデルの当てはまり具合を自由度調整済み決定係数で比べます。
消費関数の方は 0.882、投資関数の方は 0.549 なので、消費関数の方が当てはまりが良いと言えそうです。
係数の$${P}$$値を比べます。
消費関数の方は説明変数は 0.0 % であり 5%水準で有意ですが、切片は 9.9%であり 5%水準で有意ではありません。
消費関数の方は全ての係数で 0.0 % であり 5%水準で有意です。
今回の写経は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！
１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。
２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で
書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！
３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で
書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！
４．楽しい写経 ベイズ・Python等
ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀
５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で
書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。
６．データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。
７．Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。