「回帰分析から学ぶ計量経済学」をPythonで写経 ～ 第３章「式の工夫」③生産関数・対数線形・ロジスティック曲線

第３章「式の工夫」

書籍の著者　山澤成康 先生

この記事は、書籍「回帰分析から学ぶ計量経済学」第３章「式の工夫」の Python写経活動 を取り扱います。

今回は第３章の課題を実践します。
生産関数、対数線形、ロジスティック曲線に取り組みます！
では書籍を開いて回帰分析の旅に出発です🚀

はじめに

書籍「回帰分析から学ぶ計量経済学」のご紹介

このシリーズは書籍「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」（オーム社、「テキスト」と呼びます）の Python 写経です。

テキストは、2023年11月に発売され、副題「Excelで読み解く経済のしくみ」のとおり、主に Excel を用いて、計量経済学を平易に学べる素晴らしい書籍です。
テキストの「はじめに」に著者の先生が執筆の動機を書かれています。

私にとって計量経済学は高嶺の花ですが、このテキストでさまざまな回帰分析のアプローチを知ることができました。
また、書籍の Excel 処理を Python に置き換える「寄り道写経」の実践を通じて、回帰分析のお気持ちに少し近づけた感じがいたします。

回帰分析に慣れ親しむのに丁度良いレベル感と内容ですので、これはぜひともブログにしたい！と思って現在に至ります。
計量経済学の色を薄め、データ分析の色を濃いめに書いてまいります！

引用表記

この記事は、出典に記載の書籍に掲載された文章と配布データを引用し、適宜、掲載文章・配布データを改変して書いています。
【出典】
「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」
第１版第１刷、著者 山澤成康、オーム社

記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！

第３章 式の工夫

この記事は第３章の以下の節を取り扱います。

記事に用いるデータは、オーム社の書籍紹介サイトからダウンロードできる 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')

課題1 生産関数

生産関数を用いて回帰分析を実践する課題です。
テキストの生産関数の数式展開をお借りした後に、回帰分析に進みます。

１．生産関数
$${Y}$$は実質GDP、$${A}$$は全要素生産性、$${K}$$は資本ストック、$${L}$$は労働投入量です。

２．両辺を対数にして式変形

３．置き換え
以下の変数で置き換えます。
最後の数式を回帰分析に用います！

４．係数の推定値を用いた対数実質GDPの推定式

それでは回帰分析に進みます！

■ データの読み込み
テキストのデータを読み込みます。
経済指標値の出典はテキストをお読み下さい！

### データの読み込み

# CSVファイルの読み込み
df7 = pd.read_csv('./data/03_ex1.csv', index_col=0)
# データフレームの表示
print('df7.shape:', df7.shape)
display(df7.head(3))
display(df7.tail(3))

【実行結果】
1994～2022年のデータです。
欠損値がありますね・・・。

■ データの前処理
欠損値を削除し、労働投入量$${L}$$、目的変数$${y}$$、説明変数$${x}$$、トレンド変数$${T}$$を作成します。

### データの前処理

## 準備
# データフレームのコピーを取得
df7_mod = df7.copy()
# 欠損値を含む行を削除
df7_mod = df7_mod.dropna()

## 前処理
# 労働投入量Lの算出
df7_mod['労働投入量'] = df7_mod['総労働時間'] * df7_mod['就業者数']
# 目的変数yの算出
df7_mod['y'] = np.log(df7_mod['実質GDP']) - np.log(df7_mod['労働投入量'])
# 説明変数xの算出
df7_mod['x'] = np.log(df7_mod['資本ストック']) - np.log(df7_mod['労働投入量'])
# トレンド変数の作成
df7_mod['trend'] = np.arange(1, len(df7_mod)+1)

## 結果の表示
print('df7_mod.shape: ', df7_mod.shape)
display(df7_mod.head(3))
display(df7_mod.tail(3))

【実行結果】
トレンド変数は初年度が$${1}$$、以降の年度は$${+1}$$ずつ増えています。

■ 回帰分析の実行

### 回帰分析の実行

result_production = smf.ols(formula='y ~ x + trend', data=df7_mod).fit()
display(result_production.summary())

【実行結果】
係数は有意です。

対数実質GDPの推定式は次のようになります。

■ お楽しみ
実質GDP$${Y}$$の観測値と予測値を時系列で描画して締めましょう。

### 実質GDPの観測値と予測値のプロット

## 実質GDPの予測値の算出
# 回帰係数の取得
a, b, c = result_production.params
# 実質GDPの予測値の算出
gdp_pred = np.exp(a + b * np.log(df7_mod['労働投入量'])
                  + (1 - b) * np.log(df7_mod['資本ストック'])
                  + c * df7_mod['trend'])

## 描画
# 描画領域の設定
plt.figure(figsize=(10, 4))
# 観測値の散布図の描画
plt.scatter(df7_mod.index, df7_mod['実質GDP'], alpha=0.7, label='観測値')
# 予測値の散布図の描画
plt.plot(df7_mod.index, gdp_pred, color='tab:red', ls='--',
         label='予測値')
# 修飾：x軸ラベル、y軸ラベル、グリッド線、凡例表示
plt.xlabel('年度', fontsize=12)
plt.ylabel('実質GDP', fontsize=12)
plt.grid(lw=0.5)
plt.legend();

【実行結果】
まずまずの予測値ではないでしょうか。

課題2 対数線形

「いか」と「たこ」の需要の所得弾力性と価格弾力性を計算する課題です。
目的変数と説明変数を対数にして推定する対数線形を用います。
テキストの101～102ページによると、対数線形で推定する係数は「弾力性」（説明変数の伸び率に対する目的変数の伸び率）を表すそうです。
テキストの対数線形の数式をお借りします。

■ データの読み込み
テキストのデータを読み込みます。
データの出典は総務省の「家計調査」と「消費者物価指数」です。

### データの読み込み
# 購入数量：g、実収入：円、CPI：%

# CSVファイルの読み込み
df8 = pd.read_csv('./data/03_ex2.csv', index_col=0)
# 実質実収入列の追加 実収入 / 総合CPI * 100
df8['実質実収入'] = df8['実収入'] / df8['総合CPI'] * 100
# インデックスの「年」の文字を削除
df8.index = df8.index.str.replace('年', '').astype(int)
# データフレームの表示
print('df8.shape:', df8.shape)
display(df8.head(3))
display(df8.tail(3))

【実行結果】
2000～2022年のいか・たこ購入数量等のデータです。

■ データの前処理
各変数を対数にします。

### データの前処理

# 対数化
df8_log = df8.apply(np.log)
# 列名の末尾に_logを設定
df8_log = df8_log.add_suffix('_log', axis=1)
# 結果の表示
print('df8_log.shape:', df8_log.shape)
display(df8_log.head())

【実行結果】

■ いかの回帰分析
対数にした変数を用いて回帰分析を実行します。

### いかの回帰の実行
result_ika = smf.ols(
    formula='いか購入数量_log ~ 実質実収入_log + いかCPI_log', data=df8_log).fit()
display(result_ika.summary())

【実行結果】
係数は有意です。
所得の係数は$${4.8478}$$、物価の係数は$${-2.2563}$$です。

実は対数にした変数を作らなくても、statsmodels の式上で対数にする方法もあります。

### いかの回帰の実行 ※formula内で対数化する
result_ika2 = smf.ols(
    formula='I(np.log(いか購入数量)) ~ I(np.log(実質実収入)) + I(np.log(いかCPI))',
    data=df8).fit()
display(result_ika2.summary())

【実行結果】
最初の回帰分析と同じ結果になります。

いかの購入数量の観測値と予測値の時系列で可視化します。

### 観測値と予測値のプロット 対数から元の単位[g]に変換して可視化

# 描画領域の設定
plt.figure(figsize=(10, 4))
# 観測値の散布図の描画
plt.scatter(df8.index, df8['いか購入数量'], alpha=0.7, label='観測値')
# 予測値の散布図の描画
plt.plot(df8.index, np.exp(result_ika.fittedvalues), color='tab:red', ls='--',
         label='予測値')
# 修飾：x軸ラベル、y軸ラベル、グリッド線、凡例表示
plt.xlabel('年度', fontsize=12)
plt.ylabel('いか購入数量 [g]', fontsize=12)
plt.xticks(ticks=df8.index, rotation=30)
plt.grid(lw=0.5)
plt.legend();

【実行結果】
2013年までの予測値は少々外れ気味な感じがします。

■ たこの回帰分析
対数にした変数を用いて回帰分析を実行します。

### たこの回帰の実行
result_tako = smf.ols(
    formula='たこ購入数量_log ~ 実質実収入_log + たこCPI_log', data=df8_log).fit()
display(result_tako.summary())

【実行結果】
係数は有意です。
所得の係数は$${3.0715}$$、物価の係数は$${-1.6086}$$です。

たこもstatsmodels の式上で対数にする方法を試します。

### たこの回帰の実行 ※formula内で対数化する
result_tako2 = smf.ols(
    formula='I(np.log(たこ購入数量)) ~ I(np.log(実質実収入)) + I(np.log(たこCPI))',
    data=df8).fit()
display(result_tako2.summary())

【実行結果】
最初の回帰分析と同じ結果になります。

たこの購入数量の観測値と予測値の時系列で可視化します。

### 観測値と予測値のプロット 対数から元の単位[g]に変換して可視化

# 描画領域の設定
plt.figure(figsize=(10, 4))
# 観測値の散布図の描画
plt.scatter(df8.index, df8['たこ購入数量'], alpha=0.7, label='観測値')
# 予測値の散布図の描画
plt.plot(df8.index, np.exp(result_tako.fittedvalues), color='tab:red', ls='--',
         label='予測値')
# 修飾：x軸ラベル、y軸ラベル、グリッド線、凡例表示
plt.xlabel('年度', fontsize=12)
plt.ylabel('たこ購入数量 [g]', fontsize=12)
plt.xticks(ticks=df8.index, rotation=30)
plt.grid(lw=0.5)
plt.legend();

【実行結果】
急激な増減の年の予測値は外れがちですが、それ以外の年はまずまずのように感じます。

課題3 ロジスティック曲線

2人以上世帯の食洗機普及率の将来予測に取り組みます。
ロジスティック曲線を当てはめます。

■ データの読み込み
テキストのデータを読み込みます。

### データの読み込み
#普及率 %

# CSVファイルの読み込み
df9 = pd.read_csv('./data/03_ex3.csv')
# データフレームの表示
print('df9.shape:', df9.shape)
display(df9.head(3))
display(df9.tail(3))

【実行結果】
2005～2023年の普及率です。

■ ロジスティック曲線の推定
テキスト104ページのロジスティック曲線の数式をお借りします。
この数式を使って変数を変換し、$${Z}$$を目的変数にして通常の回帰分析を実行します。

### ロジスティック曲線の推定

## 回帰で使用する変数の作成
# 飽和点Sの設定
S = 100
# 回帰分析用のデータフレームの作成、列名をX,Yに変更
df9_logit = df9.copy()
df9_logit.columns = ['X', 'Y']
# 回帰分析用の目的変数Zの算出
df9_logit['Z'] = -np.log((S - df9_logit['Y']) / df9_logit['Y'])

## 回帰分析 ※ Z = 1 + X
result_logit = smf.ols(formula='Z ~ X', data=df9_logit).fit()
result_logit.summary()

【実行結果】

■ 予測の実行
2005～2050年までの普及率を予測します。
回帰モデルで$${\hat{Z}}$$を予測し、$${\hat{Y} = S / (1 + \exp(-\hat{Z}))}$$で普及率$${Y}$$の予測値を算出します。

### 予測の実行

# 予測対象年2005～2050の設定
X_pred = np.arange(2005, 2051)

# Zの推定値を算出 ※回帰モデルで予測
Z_pred = result_logit.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())

【実行結果】

■ 普及率の予測値の時系列プロット
予測結果を可視化しましょう。

### 可視化

# 描画領域の設定
plt.figure(figsize=(10, 4))
# 再生可能エネルギー比率の実績値の折れ線グラフの描画
plt.plot(df9['年'], df9['食器洗い機普及率'], '-o', label='実績値')
# 再生可能エネルギー比率の推定値の点線グラフの描画
plt.plot(X_pred, Y_pred, color='tab:red', ls='--', label='推定値')
# 修飾
plt.xlabel('年', fontsize=12)
plt.ylabel('食器洗い機普及率 [%]', fontsize=12)
plt.title('食器洗い機普及率の予測')
plt.grid(lw=0.5)
plt.legend()
plt.show()

【実行結果】
普及率の予測値は直線的になりました。

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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の教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。