「回帰分析から学ぶ計量経済学」をPythonで写経 ～ 第３章「式の工夫」①さまざまな式の形
第３章「式の工夫」
書籍の著者 山澤成康 先生
この記事は、書籍「回帰分析から学ぶ計量経済学」第３章「式の工夫」の Python写経活動 を取り扱います。
今回は回帰モデルの右辺や左辺が変形します！
原点回帰、２次関数、双曲線、対数線形、ロジスティック曲線に取り組みます！
では書籍を開いて回帰分析の旅に出発です🚀
はじめに
書籍「回帰分析から学ぶ計量経済学」のご紹介
このシリーズは書籍「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」（オーム社、「テキスト」と呼びます）の Python 写経です。
テキストは、2023年11月に発売され、副題「Excelで読み解く経済のしくみ」のとおり、主に Excel を用いて、計量経済学を平易に学べる素晴らしい書籍です。
テキストの「はじめに」に著者の先生が執筆の動機を書かれています。
社会人の統計リテラシーの向上をテーマの１つとした科研費プロジェクトの最終年度で、広く社会人に向けてわかりやすい経済分析の本を書きたかったのです。
私にとって計量経済学は高嶺の花ですが、このテキストでさまざまな回帰分析のアプローチを知ることができました。
また、書籍の Excel 処理を Python に置き換える「寄り道写経」の実践を通じて、回帰分析のお気持ちに少し近づけた感じがいたします。
回帰分析に慣れ親しむのに丁度良いレベル感と内容ですので、これはぜひともブログにしたい！と思って現在に至ります。
計量経済学の色を薄め、データ分析の色を濃いめに書いてまいります！
引用表記
この記事は、出典に記載の書籍に掲載された文章と配布データを引用し、適宜、掲載文章・配布データを改変して書いています。
【出典】
「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」
第１版第１刷、著者 山澤成康、オーム社
記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！
第３章 式の工夫
この記事は第３章の以下の節を取り扱います。
3.2 さまざまな式の形
3.3 定数項なしの推定
3.4 2次式の推定（スマイルカーブ）
3.5 トレードオフを表すグラフ（双曲線）
3.8 ロジスティック曲線
記事に用いるデータは、オーム社の書籍紹介サイトからダウンロードできる 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
# トービット、ヘーキットモデル
from py4etrics.tobit import Tobit
# 描画
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'
# ワーニング非表示
import warnings
warnings.simplefilter('ignore')
3.1 さまざまな式の形
テキストの数式とチャートを紹介する形式でお話を進めます。
■ 1次式（線形）のモデル
テキストの数式を引用し、グラフを描画します（他のモデルも同じ）。
$$
Y_i = \alpha + \beta X_i + u_i
$$
### １次式 Y = a + bX
## 設定と準備
# 係数の設定
a, b = 1, 0.5
# xの値の設定
X = np.linspace(1, 22, 1001)
# Yの算出
Y = a + b * X
## 描画
plt.plot(X, Y)
plt.title(r'$Y = \alpha + \beta X$')
plt.ylim(0, 14)
plt.grid(lw=0.5);【実行結果】
■ 原点回帰
定数項のない1次式です。
$$
Y_i = \beta X_i + u_i
$$
### 原点回帰 Y = aX
## 設定と準備
# 係数の設定
a = 0.5
# xの値の設定
X = np.linspace(1, 22, 1001)
# Yの算出
Y = a * X
## 描画
plt.plot(X, Y)
plt.title(r'$Y = \alpha X$')
plt.grid(lw=0.5);【実行結果】
■ 2次関数
説明変数が 2次式になっています。
$$
Y_i = \alpha + \beta_{1i} X_i + \beta_{2i} X^2_i
$$
（注）2次関数以降、テキストの数式には誤差項がありません。
### 2次曲線 Y = a + b1X + b2X²
## 設定と準備
# 係数の設定
a, b1, b2 = 1, 0.5, 0.3
# xの値の設定
X = np.linspace(1, 22, 1001)
# Yの算出
Y = a + b1 * X + b2 * X**2
## 描画
plt.plot(X, Y)
plt.title(r'$Y = \alpha + \beta_1 X + \beta_2 X^2$')
plt.grid(lw=0.5);【実行結果】
■ 双曲線
説明変数の逆数を使います。
$$
Y_i = \alpha + \beta (1 / X_i)
$$
### 双曲線 Y = a + b (1/X)
## 設定と準備
# 係数の設定
a, b = 1, 40
# xの値の設定
X = np.linspace(1, 22, 1001)
# Yの算出
Y = a + b * (1 / X)
## 描画
plt.plot(X, Y)
plt.title(r'$Y = \alpha + \beta\ \frac{1}{X}$')
plt.ylim(0, 45)
plt.grid(lw=0.5);【実行結果】
■ 対数線形
目的変数と説明変数を対数にします。
$$
\log(Y_i) = \alpha + \beta \log(X_i)
$$
### 対数線形 log Y = a + b log X
## 設定と準備
# 係数の設定
a, b = 2.5, 3.3
# xの値の設定
X = np.linspace(1, 22, 1001)
# Yの算出
log_Y = a + b * np.log(X)
## 描画
plt.plot(X, log_Y)
plt.title(r'$log(Y) = \alpha + \beta log(X)$')
plt.ylim(0, 14)
plt.grid(lw=0.5);【実行結果】
■ ロジスティック曲線
詳細は 3.8 節で解説されます。
$$
Y_i = \cfrac{1}{1 + e^{-X_i}}
$$
### ロジスティック曲線 Y = 1 / (1 + e^-x)
## 設定と準備
# xの値の設定
X = np.linspace(-10, 10, 1001)
# Yの算出
Y = 1 / (1 + np.exp(-X))
## 描画
plt.plot(X, Y)
plt.axhline(0.5, color='black', ls='--', lw=0.8)
plt.title(r'$Y = 1\ /\ (1 + e^{-X})$')
plt.grid(lw=0.5);【実行結果】
3.3 定数項なしの推定
回帰式に定数項が無い「原点回帰」の深堀りです。
テキストによると、予測に原点回帰を使うのは、説明変数$${X}$$と目的変数$${Y}$$の符号を一致させたい場合、とのことです。
テキストの実質GDPデータを読み込みます。
### データの読み込み
# CSVファイルの読み込み
df1 = pd.read_csv('./data/03_01_CI.csv', index_col=0)
# 四半期日付のインデックスを設定
df1.index = pd.date_range(start='1994-03-31', periods=df1.shape[0], freq='QE')
df1.index.name = '四半期'
# GDPとCIの前期比列を追加
df1['実質GDP（前期比伸び率）'] = df1['実質GDP（水準）'].pct_change() * 100
df1['CI（前期比伸び率）'] = df1['CI（水準）'].pct_change() * 100
# データフレームの表示
print('df1.shape:', df1.shape)
display(df1.head())【実行結果】
テキストにならって、目的変数を実質GDPの前期比、説明変数を景気動向指数CIの前期比とする原点回帰を行います。
■ 実質GDP前期比伸び率とCI前期比伸び率 p.95
散布図を描画して回帰直線が原点を通ることが妥当かどうかを確認します。
### GDPとCIの前期比伸び率の散布図の描画
# 散布図の描画
plt.scatter(df1['実質GDP（前期比伸び率）'], df1['CI（前期比伸び率）'],
s=70, ec='white', alpha=0.7)
# 修飾
plt.xlabel('季節調整済み実質GDP前期比伸び率 [%]')
plt.ylabel('景気動向指数CI前期比伸び率 [%]')
plt.title('GDP伸び率と景気動向指数CIの伸び率')
plt.grid(lw=0.5)
plt.show()【実行結果】
回帰直線は原点$${(0,0)}$$を通っていそうですね！
原点回帰を実行します。
デフォルトが切片なしの statsmodels.api の OLS を利用します。
ちなみに OLS は最小二乗法（Ordinary Least Squares）の略称です！
### 定数項がない場合の回帰分析 smを使った。列名に（）が含まれるとエラーになるっぽい
# 回帰分析の実行
result1 = sm.OLS(endog=df1['実質GDP（前期比伸び率）'].dropna() ,
exog=df1['CI（前期比伸び率）'].dropna()
).fit()
# 結果の表示
display(result1.summary())【実行結果】
テキストの分析結果と異なる統計量（$${F}$$値の$${P}値$$、自由度調整済み決定係数）がありますが、係数の推定値は一致しているので、気にしないことにしましょう・・・
statsmodels の推定値を用いて、テキストにならって回帰式を書きます！
$$
実質GDP前期比伸び率 = 0.3180 \times CI前期比伸び率
$$
ちなみに定数項のある回帰モデルはこちら！
### 定数項がある場合の回帰分析 smを使った。列名に（）が含まれるとエラーになるっぽい
# 回帰分析の実行
result2 = sm.OLS(
endog=df1['実質GDP（前期比伸び率）'].dropna() ,
exog=sm.add_constant(df1['CI（前期比伸び率）'].dropna())
).fit()
# 結果の表示
display(result2.summary())【実行結果】
定数項（切片）は$${0.1196}$$、係数は$${0.3156}$$になりました。
3.4 2次式の推定（スマイルカーブ）
スマイルカーブ（口角が上がった笑顔）を描きます。
数式はこちら！
$$
Y_i = \alpha + \beta X_i + \gamma X^2_i + u_i
$$
テキストのデータを読み込みます。
### データの読み込み
# CSVファイルの読み込み
df2 = pd.read_csv('./data/03_02_happy.csv')
# データフレームの表示
print('df2.shape:', df2.shape)
display(df2.head())【実行結果】
年齢 18～90 歳の幸福度（仮想データ）です。
テキストによると「多くの国で幸福度は18歳以降低下を続けるが、47歳ころに底を打ち、82歳以上で幸福度が最高になる」だそうです。
好転したい！
散布図を描画します。
### 散布図の描画
# 描画領域の設定
plt.figure(figsize=(6, 3))
# 散布図の描画
plt.scatter(df2['年齢'], df2['幸福度'])
# 修飾
plt.xlabel('年齢')
plt.ylabel('幸福度')
plt.ylim(0, 80)
plt.grid(lw=0.5)
plt.show()【実行結果】
たしかにスマイルの口が見えます。
最小二乗法で2次式の推定ができます。
statsmodels の formula API で ols に説明変数を与える際、「 I(年齢**2) 」と I() に変数の二乗を与えることで、2次の変数を表現できます。
### 回帰分析 Y = a + b₁X + b₂X²
# formulaを用いる場合、数式をI()の中にいれる。例えば年齢の二乗は I(年齢**2)
# 回帰分析の実行
result = smf.ols(formula='幸福度 ~ 年齢 + I(年齢**2)', data=df2).fit()
# 結果の表示
display(result.summary())【実行結果】
テキストの統計量等と一致しています。
■ スマイルカーブ p.97
散布図に合わせて予測値＝スマイルカーブを描画します。
### 散布図の描画
# 描画領域の設定
plt.figure(figsize=(6, 3))
# 散布図の描画
plt.scatter(df2['年齢'], df2['幸福度'], s=70, ec='white', alpha=0.7)
# 予測値の曲線の描画
plt.plot(df2['年齢'], result.fittedvalues, color='orangered', lw=3, ls='--')
# 修飾
plt.xlabel('年齢')
plt.ylabel('幸福度')
plt.ylim(0, 80)
plt.grid(lw=0.5)
plt.show()【実行結果】
右上がり気味のスマイルの口を描けました！
回帰式は次のようになります。
$$
Y_i = 83.0605 - 2.2865 X_i + 0.0237 X^2_i
$$
statsmodels の推定結果を用いて目的変数の予測を実行するコードは次のような感じになります。
### predictする場合、年齢を辞書型（キーはデータフレームの変数名）で渡せばよい
result.predict(exog=dict(年齢=df2['年齢'].values)).rename('予測値').to_frame()【実行結果】
3.5 トレードオフを表すグラフ（双曲線）
経済学のフィリップスカーブは物価上昇率と失業率のトレードオフの関係を示す曲線だそうです。
トレードオフは双曲線で示せるようです。
テキストの双曲線の数式をお借りします。
$$
Y_i = a + \cfrac{b}{X_i - c}
$$
テキストの失業率 vs 物価上昇率のデータを読み込みます。
### データの読み込み
# CSVファイルの読み込み
df3 = pd.read_csv('./data/03_03_phillips.csv')
# データフレームの表示
print('df3.shape:', df3.shape)
display(df3.head())【実行結果】
散布図を描画します。
### 散布図の描画
# 描画領域の設定
plt.figure(figsize=(6, 4))
# 散布図の描画
plt.scatter(df3['完全失業率'], df3['消費者物価指数上昇率'], s=70, ec='white',
alpha=0.7)
# 修飾
plt.xlabel('失業率 [%]')
plt.ylabel('消費者物価指数上昇率 [%]')
plt.xlim(0, 6)
plt.ylim(-5, 25)
plt.grid(lw=0.5)
plt.show()【実行結果】
下に凸の曲線がうっすらと見えるのではないでしょうか。
回帰分析に進みます。
テキストによると、1次式に変換しないと最小二乗法を適用できないため、$${X}$$軸の漸近線を決める：$${c=-2}$$とし、回帰式を変形します。
$$
Y_i = a + b \times \cfrac{1}{X_i +2}
$$
### 回帰分析 X軸の漸近線を決めないと１次式に変換できない ⇒ c=-2 とする
# formulaの場合、1/(完全失業率+2)をI()の中にいれる
# 回帰分析の実行
result = smf.ols(
formula='消費者物価指数上昇率 ~ I(1 / (完全失業率 + 2))', data=df3).fit()
# 結果の表示
display(result.summary())【実行結果】
テキストの統計量等と一致しています。
■ フィリップスカーブの推定 p.99
推定したモデルの予測値を散布図上にプロットしましょう。
### 散布図の描画
# 描画領域の設定
plt.figure(figsize=(6, 4))
# 散布図の描画
plt.scatter(df3['完全失業率'], df3['消費者物価指数上昇率'], s=70, ec='white',
alpha=0.7)
# 予測値の曲線の描画
sort_idx = df3['完全失業率'].argsort()
plt.plot(df3['完全失業率'][sort_idx], result.fittedvalues[sort_idx],
color='tab:red', lw=2, ls='--')
# 修飾
plt.xlabel('失業率 [%]')
plt.ylabel('消費者物価指数上昇率 [%]')
plt.xlim(0, 6)
plt.ylim(-5, 25)
plt.grid(lw=0.5)
plt.show()【実行結果】
緩やかな下に凸の曲線を描けました。
3.8 ロジスティック曲線
テキストのロジスティック曲線の数式をお借りします。
$$
Y_i = \cfrac{e^{X_i}}{1 + e^{X_i}} = \cfrac{1}{1 + \frac{1}{e^{X_i}}} = \cfrac{1}{1 + e^{-X_i}}
$$
■ 指数曲線とロジスティック曲線 p.103
指数関数$${Y=e^X}$$の曲線とロジスティック曲線を比較します。
### 可視化
## データ作成
# x軸の値の作成
x = np.linspace(-10, 11, 1001)
# 指数関数 y=e^x の算出
y1 = np.exp(x)
# ロジスティック関数 y = 1 / (1 + e^-x) の算出
y2 = 1 / (1 + np.exp(-x))
## 描画
# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3))
## 指数曲線の描画
# 指数曲線の折れ線グラフの描画
ax1.plot(x, y1)
# 修飾
ax1.set(xlabel='x', ylabel='y', title='指数曲線 : $y=e^x$')
ax1.grid(lw=0.5)
## ロジスティック曲線の描画
# ロジスティック曲線の折れ線グラフの描画
ax2.plot(x, y2)
# 修飾
ax2.set(xlabel='x', ylabel='y',
title=r'ロジスティック曲線 : $y=1\ /\ (1 + e^{-x})$')
ax2.grid(lw=0.5)【実行結果】
ロジスティック曲線の$${Y}$$は$${0}$$から$${1}$$の間の値をとります。
■ ロジスティック曲線の推定 p104～105
再生可能エネルギーの比率がロジスティック曲線に沿って増加する場合の予測を回帰分析を用いて行います。
式変形と最小二乗法を用いてロジスティック回帰を実現する、みたいな感じです。
しばらく式の変形を実施します。
【式変形その１】
ロジスティック曲線の式に次の変更を加えます。
・$${X_i}$$を 1次式$${\alpha + \beta X_i}$$に変更
・目的変数$${Y}$$の最大値＝飽和点を$${1}$$から$${S}$$に変更
$$
Y_i = \cfrac{S}{1 + e^{-(\alpha+\beta X_i)}} \tag{1}
$$
【式変形その２】
式1 を「えーい」と変換して$${Z_i}$$に代入します。
$$
Z_i = -\log \left( \cfrac{S - Y_i}{Y_i} \right) =\alpha+\beta X_i \tag{2}
$$
【式変形その３】
式1 の 1次式$${\alpha + \beta X_i}$$に$${Z_i}$$を代入します。
$$
Y_i = \cfrac{S}{1 + e^{-Z_i}} \tag{3}
$$
準備完了です。
テキストの再生可能エネルギー比率データを読み込みます。
### データの読み込み
# CSVファイルの読み込み
df4 = pd.read_csv('./data/03_04_renew_energy.csv')
# データフレームの表示
print('df4.shape:', df4.shape)
display(df4.head())【実行結果】
2005～2021年までの再生可能エネルギー比率のデータです。
年列が説明変数$${X}$$、再生可能エネルギー比率が目的変数$${Y}$$です。
式2 を用いて$${Z}$$を算出し、$${Z}$$に$${X}$$を回帰させます。
飽和点$${S}$$を$${100}$$（単位：％）とします。
statsmodels の formla API の ols を利用します。
### ロジスティック曲線の推定
## 回帰で使用する変数の作成
# 飽和点Sの設定
S = 100
# 回帰分析用のデータフレームの作成、列名をX,Yに変更
df4_logit = df4.copy()
df4_logit.columns = ['X', 'Y']
# 回帰分析用の目的変数Zの算出
df4_logit['Z'] = -np.log((S - df4_logit['Y']) / df4_logit['Y'])
## 回帰分析 ※ Z = 1 + X
result = smf.ols(formula='Z ~ X', data=df4_logit).fit()
result.summary()【実行結果】
テキストの統計量等と値が異なりますが、気にしないことにしましょう。
（書籍サイトでダウンロードしたExcelファイルの値とは一致しています）
回帰式は次のようになります。
$$
再生可能エネルギー比率 = \cfrac{100}{1 + e^{-(-209.1404 + 0.1022 \times 年)}}
$$
さあ、予測に進みましょう。
上記のモデルに $${X}$$（年）を$${2005}$$～$${2040}$$年まで設定して、$${Z}$$の予測値を得ます。
そして、式3 に$${Z}$$の予測値を代入して、$${Y}$$の予測値を得ます。
## 予測の実行
# 予測対象年2005～2040の設定
X_pred = np.arange(2005, 2041)
# Zの推定値を算出 ※回帰モデルで予測
Z_pred = result.predict(exog=dict(X=X_pred))
# Yの推定値を算出 ※ Y = S / (1 + e^-z)
Y_pred = S / (1 + np.exp(-Z_pred))
# 予測値の表示
display(pd.DataFrame({'年': X_pred, 'Z': Z_pred, 'Y': Y_pred}).head(10))【実行結果】
先頭 10 年分の予測値を表示しています。
Y が再生可能エネルギー比率の予測値です。
再生可能エネルギー比率の予測値を可視化しましょう。
## 可視化
# 描画領域の設定
plt.figure(figsize=(8, 3))
# 再生可能エネルギー比率の実績値の折れ線グラフの描画
plt.plot(df4['年'], df4['再生可能エネルギー比率'], '-o', label='実績値')
# 再生可能エネルギー比率の推定値の点線グラフの描画
plt.plot(X_pred, Y_pred, color='tab:red', ls='--', label='推定値')
# 修飾
plt.xlabel('年')
plt.ylabel('再生可能エネルギー比率 [%]')
plt.title('再生可能エネルギー割合の予測')
plt.grid(lw=0.5)
plt.legend()
plt.show()【実行結果】
再生可能エネルギー比率は、2021年までの実績値（青点）に沿って、緩やかに増加する、そんな曲線（赤点線）が描かれました。
今回の写経は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。