# load libraries
pacman::p_load("tidyverse", "skimr", "huxtable", "SuperLearner", "grf", "ranger", "estimatr")概要
本稿では、Arthur M. Okunが1969年に発表した論文で主張しているOkun’s Lawを基に、失業率の景気への影響のパラメーター検証を行う。分析で用いた国々は、オーストラリア、ベルギー、カナダ、デンマーク、フィンランド、フランス、ドイツ、イタリア、日本、オランダ、ノルウェー、スウェーデン、スイス、イギリス、そしてアメリカの16か国である。分析方法は、OLSとセミパラメトリックモデルの二種類を用いた。用いたデータセットは、OLSの分析に対しては、IMFのデータを用いた。データの期間が最も長くとれた国に対しては、1950年から2020年の機関でデータを取得した。セミパラメトリックモデルは、OECDから取得したデータセットで分析をし、1985年から2011年のデータを使って検証した。Okun’sLawは発表から半世紀余りたっているが、多くの国において、過去の研究(e.g. Goto, 2020; Jim, 2000)が示すように、失業率と景気の間に負の関係性があることが知られている。本稿では分析にR言語を用いて、Okun’s Lawのパラメーターの推定を行った。
Okun(1962)が発表した論文は、失業率とビジネスサイクルにおける変動について触れており、そのパラメーターは3パーセントであることが示されていた。このパラメータが3パーセントであるというのは、景気が1パーセント変動した場合、それに応じて失業率が3パーセント負の方向に変動するということである。(JIM, 2000) 失業率と景気の負の関係性は、rule of thumとして経済界において広く受け入れられており、Mankiw(2012)などの世界的に知名度の高い経済学の教科書の中でも、実証された規則性だと説明されている。そして、Okun’s Lawは、NAIRUなどの理論的基礎となっている。
しかし、私たちは、失業率と景気の度合いを強調しすぎているのかもしれない。Josep E(2020)が主張するように、GDPそれ自体の価値を見直そうとする動きもある中、筆者は、GDPが、失業率を考える際において、本当に重要性を持っているのか疑問に考えた。本稿では失業率が景気に及ぼす程度を、主にOLSと機械学習を用いてそのパラメータを計測するという形で論を進めていく。本稿は、2章データと分析方法、が3章がそれぞれの分析結果、4章が結果に対する考察、そして5章が結論となっている。本稿の最後にAPPENDIXで式の展開の説明を加えているため、必要があれば参照してほしい。
本稿では以下の方法でOkun’s Lawのパラメータ推定を行う。
本稿では2つのデータセットを用意した。一つ目のデータセットは、Jim(2000)で使われていたOECE16か国をもとにデータセットを作った。(オーストラリア、ベルギー、カナダ、デンマーク、フィンランド、フランス、ドイツ、イタリア、日本、オランダ、ノルウェー、スウェーデン、スイス、イギリス、そしてアメリカ)。後述するThe Difference modelの方法では、IMFで取得可能であった1950年から2020年の期間のデータを使用した。使用した変数は、失業率と実質GDPである。セミパラメトリックモデルでパラメータを推定する際に使用したデータセットは、先のThe Difference modelで用いたデータセットに加え、OECDのデータベースから、資本投入量、労働時間と労働生産性のデータである。そのため、先のデータセットと異なり、データが利用可能であった期間は一番長い国でも、1985年から2011年の27年間であった。欠損値が生じた際には、各国ごとに欠損値の修正を行い、その欠損値を含む列ごとデータセットから削除した。以下がデータフレームの全体の情報である。
# set the path
path <- "~/R/Musashi_University/2021_second_semester/final-thesis/Data/df_imf"
# import data as df
data <- read_rds(path) %>%
na.omit()
# show
DT::datatable(data,
rownames = FALSE,
extensions = 'Buttons',
options = list(autoWidth = TRUE,
pageLength = 5,
dom = 'Bfrtip',
buttons = list("csv"),
scrollX = TRUE,
scrollCollapse = TRUE),
class = 'cell-border stripe',
caption = "source: ")data2 <- read_rds("~/R/Musashi_University/2021_second_semester/final-thesis/Data/all_oecds.rds")
# show
DT::datatable(data2,
rownames = FALSE,
extensions = 'Buttons',
options = list(autoWidth = TRUE,
pageLength = 5,
dom = 'Bfrtip',
buttons = list("csv"),
scrollX = TRUE,
scrollCollapse = TRUE),
class = 'cell-border stripe',
caption = "source: ")GDPが失業率に影響する度合いを測るために、本稿では2つの方法を用いる。最初の一つは伝統的経済学的方法である、変数と関数系の選択を行ってから、それを基に回帰分析を回す方法である。具体的には、Jim(2000)で提案されていたDifference modelを用いる。Jim(2000)で議論されていたThe Gap Modelは、HP filter、BN filterなどのdetrendingの方法がが自身のリサーチ能力を超えるため、本稿では省略する。それぞれの式の導出は以下のとおりである。
この方法では、経済の総生産量(y_t)と失業数(u_t)の関係性は以下のように定義される。
\[ \Delta y_t= β_0- β_1 \Delta u_t+ ϵ_t, \qquad t=1,…,T, \qquad (1) \]
この式では、\(\Delta\)は、変化量を表しており、\(\Delta u_t\)と\(\Delta y_t\)はそれぞれGDPと失業率の変化率を示している。ϵは、攪乱項を示している。パラメータβ_0は、この式の切片を表しており、β_1は、この式の傾きを表している。 しかし、The Difference methodはトレンドを考慮していないため、パラメータにバイアスがかかってしまう。そして、この式は、二変数が線形であるという前提のもとにGDPの成長率の失業率に対する影響を計測しているため、モデルが誤って定義されている可能性が高い。そこで、二つ目に用いる方法は、機械学習を用いて、パラメータの推定を行う。
機械学習を使い、Robinsson推定によって失業率が景気に及ぼす係数を推定した。式は、以下のように設定した。
\[ y_t = τDt + f(x) + \epsilon_t,\qquad t=1,…,T \qquad (2) \]
先ほどと同じく、\(y_t\)はその年のGDP成長率を表しているが、\(D_t\)は、その年の失業率の変化率を表している。今回の推定では、\(\tau\)が目的関数であり漸近正規性を満たすと考える。\(f(x)\)はコントロール変数である。OLSのような伝統的な計量経済学の方法と異なり、\(f(x)\)にはどのような関数形が来るのかは完全にブラックボックスである。\(f(x)\)中には、Okun(1962)で提案されている変数、資本投入量、労働時間と労働生産性をコントロール変数として用いた。これらをコントロール変数として選んだ理由を以下に記述する。
まず、経済が不況であるとき、一般的に企業は利益が減少する。生産要素が資本と労働に限定され、他の条件一定の場合、企業が経済不況に対処する方法は三つある。一つ目は、資本を下げて対応する方法である。そのため、資本の変動を考慮するため、一人当たりの資本の成長率をモデルに組み入れた。
次に、不景気の際、企業は労働力を変動させて対応できる。労働力の変動は、主に雇用数を調整する方法と、各労働者の労働時間を減少させる二つの方法がある。雇用の変動は、失業率としてモデルに組入れられているが、労働時間の変動は、考慮されていないため、年間の労働時間の成長率をコントロール変数として取り入れた。
最後に、企業は不景気の際、資本も労働も変動させないという対処方法を選択できる。この場合、企業の生産性が生産量/資本の投入量で計算されるとき、生産量のみの減少となり、生産性が悪化する。不景気化における企業の生産性の悪化をモデルに取り入れるために、MFP(Multifactor Productivity)をモデルに取り入れた。
以上、不景気の際に企業が取れる行動を考慮するために、資本投入量、労働時間そして労働生産性をコントロール変数としてモデルに取り入れた。
Robinsson推定のステップは主に4つある。まず、\(y_t\)と\(D_t\)についての予測を機械学習を用いて行う。その後、交差検証を用いて、モデルの性能評価を行う。そして、\(y_i - f_y(X_i)\)と\(D-f_D(x_i)\)の単回帰をおこない、\(\hat tau\)の推定を行う。このままでは、\(\hat tau\)は、各年におけるパラメータであるため、最後に データセットから予測された\(\hat tau\)の平均と分散から、95%信頼区間を計算して、OCEDq6か国、各国の景気における失業率のパラメータを予測する。以下に各ステップについて詳細に記述する。
まず、\(y_t\)と\(D_t\)の予測には、OLS、ランダムフォレスト、LASSOを組み合わせて予測を行うstackingという方法を用いて、パラメータの予測値を求めた。OLSは、最小二乗法のことをさしており、The Difference modelと変わらない。しかし、決定木をもちいて学習を行うRandom Forestと、OLSに制約条件を加えたモデルであるLASSOと組み合わせることでパラメータの予測向上を図った。主に、パッケージSuperLeanerとパッケージgrfを使い、モデルの学習を行わせた。各モデルがどのようなMSEを最小化しようとしているのかは、以下のとおりである。
\[ Y_i = \beta_0 + \beta_1X_2 + ... + \beta_n X_n + \epsilon_i \]
\[ \sum_{i = 1}^{i = n}{u^2} \\ = \sum_{i = 1}^{i = n}{(Y - \hat Y)^2} \\ = \sum_{i = 1}^{i = n}{(Y - \hat \beta_0 - \hat \beta_1 X_2 - ... - \hat \beta_n X_n)^2} \\ \]
ランダムフォレストは、説明変数に離散変数を用いる。今回扱う変巣は、すべて連続変数だが、GDPの成長率が2パーセント以上と以下に区切るなどして、連続変数を離散変数として使用するため、連続変数を説明変数として用いても問題がない。
\[ Y_i = E[Y_i|X_i = x] + \epsilon_i \]
\[ \sum_{i = 1}^{i = n}{u^2} \\ = \sum_{i = 1}^{i = n}{(Y - f_Y(X_i))^2} \\ = \sum_{i = 1}^{i = n}{(Y - E[Y_i|X_i = x])^2} \\ \]
LASSOは、OLSの過学習を防ぐために、\(\beta\)がとれる値に制約条件を加え、MSEを最小化したモデルである。
\[ Y_i = \beta_0 + \beta_1X_2 + ... + \beta_L X_L + \epsilon_i \]
\[ \sum_{i = 1}^{i = n}{u^2} \qquad s.t. \sum_{l = 1}^{L}|B_l| \leq A \\ = \sum_{i = 1}^{i = n}{(Y - f_Y(X_i))^2} \qquad s.t. \sum_{l = 1}^{L}|B_l| \leq A \\ = \sum_{i = 1}^{i = n}{(Y - \hat \beta_0 - \hat \beta_1 X_2 - ... - \hat \beta_l X_l)^2} \qquad s.t. \sum_{l = 1}^{L}|B_l| \leq A \\ \]
これらの3つのアルゴリズムをstackingにより組み合わせ、最善の\(f_Y(X_i)\)と\(f_D(X_i)\)の報告をさせる。
次に、交差検証を行い、MSEを最小化するパラメータを設定する。今回のデータの分割数は10に設定した。
その後、\(y_i - f_y(X_i)\)と\(D-f_D(x_i)\)の単回帰をおこない、今回の目的関数である\(tau\)の推定を行う。単回帰の式は以下のようになっている。
\[ Y_i - f_Y(X_i) = \tau (D_i - f_D(X_i)) + \epsilon \]
先ほどのOLSを使うThe Difference modeltと比べて、この機械学習を用いるセミパラメトリックモデルには利点がある。一つは、コントロール変数の関数系を線形としないことでより柔軟に想定できる。そして、パラメータ\(\hat \tau\)を部分線形として求めるため、目的関数と被説明変数の関係性が明示的であることがあげられる。部分線形用いず、すべての説明変数に対してロビンソン推定のみで求める場合、各説明変数と被説明変数の関係性が複雑になるため、一見しただけでは理解できない。加えて、OLSと違い、セミパラメトリックモデル正しいモデルでなくとも、\(\hat \tau\)が漸近正規性を満たすため、社会科学での実現可能性が高い。以上の利点から、今回は機械学習を用いてパラメータの推定を行うセミパラメトリックモデルを用いた。
最後に、データが取得できた期間を通しての失業率の景気に対するパラメータを推定するため、\(\hat \tau\)の単純平均を求めた。以下がAIPWの式である。
\[ E[Y_i|D_1, x] + \frac{D_1 x (Y_i - E[Y_i | D_1])}{E[D_1|x]} + ... + E[Y_i|D_n, x] + \frac{D_n x (Y_i - E[Y_i | D_n])}{E[D_n|x]} \]
# create ordered data frame
df_order <- read_rds("~/R/Musashi_University/2021_second_semester/final-thesis/Data/unemployment_benefits.rds")
# modeling of the difference model
df_lm <- data %>%
group_by(country) %>%
mutate(gdp_rate = c(NA, diff(log(gdp), lag = 1))) %>% # log difference
nest() %>%
left_join(df_order) %>% # add the `ranking` variable
mutate(model = map(data, ~lm_robust(gdp_rate ~ uem_rate, data = .)))
# define the function for plotting
errorbar_lm <- function(data, rank = "", title = "", subtitle = "", caption = "") {
data$model <- data$model %>% set_names(c(data$country %>% unique()))
data$summary <- lapply(data$model, summary)
# create vars for plotting data
coefficient <- c()
p_value <- c()
sd_error <- c()
for (i in 1:length(data$country)) {
coefficient[i] <- data$summary[i][[1]]$coefficients[2, 1]
p_value[i] <- data$summary[i][[1]]$coefficients[2, 4]
sd_error[i] <- data$summary[i][[1]]$coefficients[2, 2]
}
# create a data frame for plotting error bar
df_graphics <- tibble(
country = data$country,
coefficient = coefficient,
p_value = p_value,
sd_error = sd_error,
ci_lower = coefficient - 1.96*sd_error,
ci_upper = coefficient + 1.96*sd_error,
p_logical = ifelse(p_value > 0.05, "p>0.05", "p<=0.05"),
value = select(data %>% ungroup, {{rank}}) %>% unname() %>% unlist()
# value = data$une_benefit
)
# plot with ggplot2
ggplot(df_graphics,
aes(value,
coefficient,
colour = p_logical)) +
geom_hline(yintercept = 0, alpha = 0.25, linetype = "dashed") +
ggrepel::geom_label_repel(aes(label = country), check_overlap = FALSE) +
geom_point() +
geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper)) +
scale_color_manual(values = c( "#33ccff", "#ff9980")) +
labs(title = title,
subtitle = str_wrap(subtitle, 80),
x = "Country",
y = "Coefficient",
caption = caption) +
# facet_wrap(~p_logical, scales = "free_x") +
theme_minimal() +
theme(legend.position = "none")
}
# plot error bars
errorbar_lm(df_lm,
rank = une_benefit,
title = "Beta1 Accross OECD countries with the difference model",
subtitle = "Ordering x-axis by the unemployment benefits",
caption = "Graph1") 上記のGraph1は、縦軸が国、横軸が\(\beta_1\)の値を95%信頼区間で示し、その国々をOECDの2019年時点での所得に対する失業手当の割合をもとに並び変えたものである。青く示されている国々は、95%信頼区間において、\(\beta_1\)の値が0になるという帰無仮説が棄却される国々である。デンマーク、フランス、ノルウェー、スウェーデンをのぞく12か国では、p値が、95%信頼区間において0.05以下になった。一方、極僅かであるが、赤く示されている国では、p値が0.05以上になる、つまり、\(\beta_1\)が、0になるという帰無仮説を棄却することができなかった国々である。各国についての詳しいデータは以下のRegression tableを参考にしてほしい。
# show regression table
first_regtable <- df_lm$model %>%
set_names(c(data$country %>% unique())) %>%
head(8)
second_regtable <- df_lm$model %>%
set_names(c(data$country %>% unique())) %>%
tail(8)
# show
huxreg(first_regtable)| Austria | Australia | Belgium | Canada | Switzerland | Germany | Denmark | Finland | |
|---|---|---|---|---|---|---|---|---|
| (Intercept) | 0.017 *** | 0.030 *** | 0.014 ** | 0.028 *** | 0.017 *** | 0.011 ** | 0.020 *** | 0.022 *** |
| (0.003) | (0.002) | (0.004) | (0.003) | (0.003) | (0.004) | (0.002) | (0.003) | |
| uem_rate | -0.001 * | -0.001 *** | -0.001 * | -0.001 *** | -0.000 | -0.001 ** | -0.000 | -0.001 *** |
| (0.001) | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | |
| N | 25 | 40 | 25 | 48 | 40 | 29 | 46 | 40 |
| R2 | 0.441 | 0.605 | 0.155 | 0.493 | 0.018 | 0.281 | 0.404 | 0.548 |
| *** p < 0.001; ** p < 0.01; * p < 0.05. | ||||||||
| France | United Kingdom | Italy | Japan | Netherlands | Norway | Sweden | United States | |
|---|---|---|---|---|---|---|---|---|
| (Intercept) | 0.020 *** | 0.017 *** | 0.002 | 0.007 | 0.016 *** | 0.027 *** | 0.024 *** | 0.032 *** |
| (0.003) | (0.003) | (0.005) | (0.003) | (0.003) | (0.002) | (0.003) | (0.002) | |
| uem_rate | -0.000 | -0.001 ** | -0.001 * | -0.001 * | -0.001 *** | -0.000 | -0.001 *** | -0.001 *** |
| (0.000) | (0.000) | (0.001) | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | |
| N | 49 | 24 | 25 | 26 | 25 | 47 | 27 | 69 |
| R2 | 0.041 | 0.657 | 0.162 | 0.367 | 0.557 | 0.071 | 0.476 | 0.619 |
| *** p < 0.001; ** p < 0.01; * p < 0.05. | ||||||||
# define the function for the Robinsson Estimate
machine_learn_estimate <- function(country, v = 2) {
data <- data2 %>%
filter(country == {{ country }}) %>%
mutate(gdp_rate = c(NA, diff(log(gdp), lag = 1))) %>% # the log difference
na.omit()
# parameter estiamation with ML -----------------------------------------------------
# set variables
Y <- data$gdp_rate # outcome var
D <- data$uem_rate # objective var
X <- select(data, hrsg, cap, mfp) # other variables
# how many of divisions?
v <- v
# used libraries
SL.libraries <- c("SL.lm", # OLS
"SL.ranger", # Random Forest
"SL.glmnet") # LASSO
# Fitting
fits <- list()
fits$fit.Y <- CV.SuperLearner(Y = Y,
X = X,
SL.library = SL.libraries,
V = v # how many of divisions?
)
fits$fit.D <- CV.SuperLearner(Y = D,
X = X,
SL.library = SL.libraries,
V = v # how many of divisions?
)
# estimate tau ------------------------------------------------------------
# rename the dataset
data_result <- data2
# 予測誤差
Y.oht <- fits$fit.Y$SL.predict
D.oht <- fits$fit.D$SL.predict
X.matrix <- model.matrix(~ 0 + .,
data = X)
# prediction for tau of D
fit <- causal_forest(Y = Y, # outcome var: working hour
W = D, # explanatory var: gender
X = X.matrix, # control var: age, exp
Y.hat = Y.oht, # computated value of outcome var
W.hat = D.oht) # computated value of explanatory var
pred <- predict(fit,
estimate.variance = T) # show the estimate variance
# AIPWの推定 (データ全体としての予測の正確性についての検証)
ate <<- average_treatment_effect(fit,
target.sample = "overlap")
df_ate <<- tibble(
country = {{ country }},
tau = ate[[1]],
ci.low = ate[[1]] - 1.96*ate[[2]],
ci.high = ate[[1]] + 1.96*ate[[2]]
)
}
# all OCED countries in the sample dataset
oecd16 <- data$country %>% unique()
# create empty data set
data_graphics <- data_frame(country = c("example"),
tau = c(1),
ci.low = c(1),
ci.high = c(1))
for(country in oecd16) {
machine_learn_estimate(country, v = 10)
data_graphics <- data_graphics %>%
add_row(country = df_ate$country, # country
tau = df_ate$tau, # estimated tau
ci.low = df_ate$ci.low, # lower confidence interval
ci.high = df_ate$ci.high) # higher confidence interval
}
# eliminate the first column
data_graphics <- data_graphics %>%
filter(country != "example") %>%
mutate(p_logical = if_else(ci.high > 0, "p>0.05", "p<=0.05"))
# visualisation
ggplot(data_graphics, aes(country,
tau,
ymin = ci.low,
ymax = ci.high,
color = p_logical)) +
geom_hline(yintercept = 0, alpha = 0.25, linetype = "dashed") +
geom_point() +
geom_errorbar() +
scale_color_manual("", values = c("#ff9980", "#33ccff")) +
labs(title = "tau estimated by Semiparametric regression for 16 OECD countries",
x = "Country",
y = "tau",
caption = "Graph2") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90))Graph2は、セミパラメトリックモデルを用いて、景気に対する失業率のパラメータを図示化したものである。\(\hat \tau\)は、すべての国において95%信頼区間において、予測値がゼロになるという帰無仮説を棄却できないという結果になった。グラフの中の点は各国の予測値であり、予測された値は、すべての国においてゼロ以下になった。これは、Okun’s Lawで提案されていた失業率と景気の間の負の関係性と一致する。
The Difference modelと機械学習を通じて得られたものをそれぞれ記述していく。まず、The Difference modelに関しては、国ごとのパラメータの違いについて言及したい。Graph1は2019年時点での所得に対する失業手当の割合で左から少ない順位で並び変えたのものである。Graph1から観察されることは、失業手当がOkun’s Lawのパラメータの国ごとの異質性との関係性が低いと考えられる。Okun’s Lawが統計的に優位でない国は、ビジネスサイクルによって起きる失業への対処の仕方によって生じるのかもしれない。筆者のリサーチの能力不足で統計的にこの並び順に意味があるのかは議論できないため、今後の課題である。また、そもそもデータが年次であったため、\(\hat \beta_1\)がうまく収束しなかったのかもしれない。
セミパラメトリックモデルでは、すべての国において、\(\hat \tau\)の値の95%信頼区間にゼロが含まれる結果となってしまった。原因として考えられることとしては、サンプル数の少なさである。今回最もデータが取得可能であった国でも年次で、26年分しか取得できず、このデータ数では\(\hat \tau\)に漸近正規性があったとしても、\(\hat \tau\)が真の値に収束する確率が低いと考えられる。
今回は、OLSのモデルについては、失業率とビジネスサイクルのみを考慮した。セミパラメトリックモデルにおいては、資本と労働力もモデルに入れたが、モデルの推定にあたってテクノロジーの発展は考慮されていない。将来的、もしくは、今現在起こっていることとしてテクノロージーの発展による構造的失業があり、それを考慮することでまたOkun’s Lawのパラメータが変動していくかもしれない。
本稿では、Okun’s Lawを基に、失業率とビジネスサイクルの関係性について、OLSとセミパラメトリックモデルを用いたパラメータ推定行った。The Difference modelからは、Okun’s Lawは多くの国において、いまだに景気と失業率の間の負の関係性を説明しているといえる。しかし、今回は、時間の移り変わりによるOKun’s Lawのパラメータの推移や、The Gap modelや、時系列データがもたらす自己相関、時系列データのトレンド処理など、取り組めなかったことも多い。それらの要素を考慮することで、Okun’s Lawのパラメータが変わっていくのかもしれない。今回取り組めなかったことは、今後の研究課題としたい。
Arthur M. Okun, (1962). Potential GNP Its measurement and significance
Joseph e., Stglitz, (2020). Measuring What Matters, SCIENTIFIC AMERICAN
Mankiw, N. G., (2012). Macroeconomics, Eighth Edition, Worth PublishersGDPを対数差分をとると成長率に近似できる導出は以下のとおりである。
\(\Delta y_t\)をt期におけるGDP成長率だと仮定したとき、
\[ \Delta y_t = \frac{Y_t-Y_{t-1}}{Y_{t-1}} \qquad (1) \]
\(log(1+ \Delta y_t)\)をテイラー展開すると、
\[ \Delta y_t \approx log(1+g_y) \qquad (2) \]
に近似される。
(1)を(2)に代入すると、
\[ \Delta y_t \approx log(1+g_y) \\ = log(1+\frac{Y_t-Y_{t-1}}{Y_{t-1}}) \\ = log(\frac{Y_t}{Y_{t-1}}) \\ = logY_t - log_{t-1} \]
よって、成長率(\(g_y\))は、その期と前期のGDPの対数差分によって計算される。