(Bing Image Creator で生成: プロンプト: Nature landscape, long view of sand desert, close up of yellow and red cactus flowers, photo.)
の続きです。今回は Log-LogモデルでのSerail Correlationの有無をチェックしたいと思います。
はじめは、Static Model です。log(tpx) ~ beta0 + beta1 * log(Lg5) + u というモデルです。tpx: TOPIX(東証株価指数), Lg5: 法人税収入【億円】です。
log(tpx) = -7.03 + 1.53 * log(Lg5) + u と推定されました。
残差を求めます。
u = beta0 + beta1 * u_lag + z のモデルを推定します。z は残差です。
u = 0.00 + 0.82 * u_lag + z と推定されました。0.82 という係数のp値はとても小さいので、u と u_lag は相関があることがわかります。つまり、static_log_resid のモデルは Serial Correlation がある、ということですね。
u と u_lag の散布図を描いてみます。
散布図でも確かに u と u_lag は相関していることがわかります。
Serial Correlation があったときは、lm() 関数で作成したモデルを summary() 関数で算出した p値 は誤りなので、Serial Correlation-Robust な p値を求めたいと思います。
Using R, Python and Julia for Introductory Econometrics
にある方法でやってみます。
まず、lmtest, sandwich の2つのパッケージを読み込みます。
そうしたら、coeftest() 関数を使います。このときに vcovHAC というのを使います。
p値 は非常に小さな値ですね。Serail Correlation を考慮しない場合は、
です。Std. Error, t value を比較すると、Serial Correlation Robust Inference のほうが値がだいぶ小さくなっていることがわかります。
信頼区間も求めておきましょう。
Serial Correlation Robust な 95%信頼区間は、1.37 ~ 1.69
普通の 95%信頼区間は、1.48 ~ 1.58 です。
ここまでの p値, 信頼区間は 理論ベース(算出方法の公式ベース)のものでした。
ここからは、シミュレーションベースの信頼区間を求めたいと思います。
まず、infer パッケージを読み込みます。
https://infer.netlify.app/articles/observed_stat_examples#two-numerical-vars---slr-1
のサイトを参考にして作業をしていきます。
df に log(tpx), log(Lg5) を追加しておきます。
次は、l_Lg5 の係数を求めます。
1.53, lm()関数の結果と同じです。
次は、ブートストラップで1000個の係数を算出します。
この1000個の推計値から信頼区間を求めます。
1.48 ~ 1.58 です。Serial Correlation を考慮してない、通常の信頼区間と同じですね。
今回は以上です。
最初から読むには、
です。
今回のコードは以下になります。
# log-logでやってみる
#
# モデルの推定
# log-log の Static Model
static_model_log <- lm(log(tpx) ~ log(Lg5), data = df)
summary(static_model_log)
#
# static_model_log の残差
u <- resid(static_model_log)
static_log_resid <- tibble(
u = u,
u_lag = c(NA, u[1:(length(u) - 1)])
)
static_log_resid
#
# u ~ u_lag + z の推定
lm(u ~ u_lag, data = static_log_resid) |>
summary()
#
# u と u_lag の散布図
ggplot(static_log_resid, aes(x = u, y = u_lag)) +
geom_point()
#
# lmtest, sandwich の読み込み
library(lmtest)
library(sandwich)
#
# Serial Correlation Robust Inference
coeftest(static_model_log, vcovHAC)
#
# 通常の Inference
coeftest(static_model_log)
#
# 95%信頼区間
confint <- tibble(
lower_robust = 1.52912 - 1.96 * 0.08301,
upper_lobust = 1.52912 + 1.96 * 0.08301,
lower_normal = 1.52912 - 1.96 * 0.026009,
upper_normal = 1.52912 + 1.96 * 0.026009
)
confint
#
# infer パッケージを読み込む
library(infer)
#
# dfに log(tpx), log(Lg5)を追加する
df <- df |>
mutate(l_tpx = log(tpx),
l_Lg5 = log(Lg5))
#
# l_Lg5 の係数を求める
slope_hat <- df |>
specify(l_tpx ~ l_Lg5) |>
calculate(stat = "slope")
slope_hat
#
# bootstrap
set.seed(360)
boot_dist <- df |>
specify(l_tpx ~ l_Lg5) |>
generate(reps = 1000, type = "bootstrap") |>
calculate(stat = "slope")
boot_dist
#
# 信頼区間を求める
percentile_ci <- get_ci(boot_dist)
percentile_ci
#
