先日、『gretlで計量経済分析』という本に載っている二段階最小二乗法の演習をRでやってみるエントリ（リンク）を書いたのですが、同じ本の次の章は「マクロ計量モデル入門」となっていて、1980年から2009年までの日本経済のデータを用い、5本の構造方程式と2本の定義式、7つの内生変数と9つの外生変数からなる超小型のマクロ計量モデルをつくるというものでした。
Rで連立方程式を解く方法の確認がてら、Rで実行してみます。
最近、研究室の学生にRの使い方をイチから教えてるのですが、これをさらに単純化して練習問題に使えるかも……？
データは本のサポートページからダウンロードできるものですが、「model.gdt」というファイルをgretlで読み込んでcsvで書き出したものを、Rに読み込ませて使いました。
準備
データは↓のようなイメージです（画像は一部です）。
library(lmtest)
library(car)
library(nleqslv)
d_macro <- read.csv('data/csv/model.csv', header=T)
colnames(d_macro)[1] <- "YEAR"
元データは1980年から2009年までの日本のマクロ経済のデータですが、最大で「2期前」を取るラグ変数があることから、教科書では全て「1982年以降」のデータでパラメータを推定しているので、ここで1、2行目はカットしておきます。
d_macro <- d_macro[3:30,]
構造方程式（推計式）群のパラメータ推定
この章では、著者によると「内生性が顕著でなく、結果があまり違わないので、全て2SLSではなくOLSで推定する」とのことなので、それに合わせます。二段階最小二乗法にするなら先日のエントリの方法になりますが、そんなに手間は変わらないです。
内生変数になる「民間最終消費支出」「民間企業設備投資」「輸入」「民間可処分所得」「利子率」を推計する回帰モデルを作っていきます。
モデルの説明もコメントで書いておきますが、たとえば民間最終消費支出は、
est1 <- lm(data=d_macro, formula=CP~YD_1 + WH_1)
こういうふうになります。
推定はもうこれでできているんですが、あとで連立方程式を解く処理のところで使いやすいように、predict()をラップして「説明変数を与えると、予測値（理論値）を吐き出す」ようにしておきます。
CP_fn <- function(YD_1, WH_1){
YD_1 <- as.numeric(YD_1); WH_1 <- as.numeric(WH_1)
prd <- predict(est1,data.frame(YD_1=YD_1, WH_1=WH_1))
return(prd)
}
as.numeric()しているのは、事情があります。後でここには、外生変数をまとめたデータフレームから特定の行を抜き出したあとに列名称で値を探して入れることにしたのですが、データフレームの列名称を要素名として持ったままの数字を引っ張ってくると、ここでdata.frame()内での列名称指定が無視されてしまいました。
なので、ごちゃごちゃするけどas.numeric()を挟むことにしました。
なお、predict()関数には、newdataという引数にデータフレーム型でデータを与える必要があります。
他の変数についても同じようにやっていきます。
est2 <- lm(data=d_macro, formula=IP ~ dGDP12 + INT_1 + IP_1)
IP_fn <- function(GDP_1, GDP_2, INT_1, IP_1){
GDP_1 <- as.numeric(GDP_1); GDP_2 <- as.numeric(GDP_2)
INT_1 <- as.numeric(INT_1); IP_1 <- as.numeric(IP_1)
prd <- predict(est2,data.frame(dGDP12=(GDP_1-GDP_2), INT_1=INT_1, IP_1=IP_1))
return(prd)
}
est3 <- lm(data=d_macro, formula=IM ~ IM_1 + GDP)
IM_fn <- function(IM_1, GDP){
IM_1 <- as.numeric(IM_1); GDP <- as.numeric(GDP)
prd <- predict(est3,data.frame(IM_1=IM_1, GDP=GDP))
return(prd)
}
GDPTAX <- d_macro$GDP - d_macro$TAX
est4 <- lm(data=d_macro, formula=YD ~ GDPTAX)
YD_fn <- function(GDP, TAX){
GDP <- as.numeric(GDP); TAX <- as.numeric(TAX)
prd <- predict(est4, data.frame(GDPTAX=(GDP-TAX)))
return(prd)
}
est5 <- lm(data=d_macro, formula=INT ~ l_GDP + l_M2CDN)
INT_fn <- function(GDP, M2CDN){
GDP <- as.numeric(GDP); M2CDN <- as.numeric(M2CDN)
prd <- predict(est5, data.frame(l_GDP=log(GDP), l_M2CDN=log(M2CDN)))
return(prd)
}
ところで、5番目の利子率の式は、教科書とはやり方を変えています。
GDPも内生変数なので、ln(GDP)で説明されるINTはGDPと同時に決まるという構造ですが、gretlでは連立方程式を解かせるときに、「内生変数として計算されたGDPのlogを取って、別の内生変数の説明変数にする」という処理が書けない（？）らしいので、教科書ではちょっと変なやり方で代用していました。
どういうやり方かというと、推定用元データ上は、ln(GDP)はエクセルか何かで予め計算してあって、元データのテーブルに一つの変数として入っているのですが、その上でまず、
ln(GDP)〜GDP
という回帰分析を行っています。（つまり、対数を取る処理を線形モデルで代用している…）
この回帰分析の結果を使って、GDPをln(GDP)ではない何かに変形したものを、疑似ln(GDP)としてINTの説明変数に入れているわけです。
しかしここではそんなことはせずに、普通に、GDPのlogをとってそれをINTの予測に用いるようにしておきます。
なので、教科書とは、最終の結果が少しずれます。（ここを教科書どおりにすると最終結果がピッタリ一致することは確認しました。）
lmの結果をまとめてCSV出力する
べつにやる必要はないのですが、上の回帰分析の結果をまとめて表にして出力する関数も書いておきます。
save_lm_results <- function(model, model_name='none'){
sum <- summary(model)
coef <- sum$coefficients
R2 <- c(sum$r.squared, NA, NA, NA)
Adj_R2 <- c(sum$adj.r.squared, NA, NA, NA)
DW_test <- c(dwtest(model)$statistic, NA, NA, dwtest(model)$p.value)
BG_test <- c(bgtest(model)$statistic, NA, NA, bgtest(model)$p.value)
results <- rbind(coef,R2, Adj_R2, DW_test, BG_test)
if(length(model$coefficients)-1 >= 2){
VIF <- c(NA, vif(model), NA, NA, NA, NA)
} else{
VIF <- c(NA, NA, NA, NA, NA, NA)
}
results <- cbind(results, VIF)
if(model_name=='none'){
output_name <- as.character(paste(model$call$formula[2], model$call$formula[1], model$call$formula[3], sep=''))
} else{
output_name <- model_name
}
write.csv(results, paste('results/lm_result_', output_name, '.csv', sep=""), row.names=TRUE, quote=FALSE)
}
lapply(X=list(est1, est2, est3, est4, est5), FUN=save_lm_results)
恒等式（定義式）群の定義
マクロモデル中の定義式（パラメータを推定するのではなく変数間の定義上の関係を記述するもの）を記述しておきます。
GDP_fn <- function(CP, IP, CG, IG, IH, ST, EX, IM){
def <- CP + IP + CG + IG + IH + ST + EX - IM
return(def)
}
GDPN_fn <- function(GDP, GDPDEF){
def <- GDP*GDPDEF/100
return(def)
}
外生変数テーブルを作る
さてここからは、推定されたパラメータをつかってシミュレーションを行っていきます。
外生変数をどうやって与えようかなと思ったのですが、なんとなく、データフレームに整理しておいて、そのうち1行を選択して与えると、全ての外生変数を読み込んで内生変数を導いてくれるような組み立てがわかりやすいかなと思い、そうしてみました。
まず、パラメータの推定に使ったデータから、外生変数として使用する変数群を抜き出してまとめます。
Exog_df <- d_macro[,c(
'YEAR',
'YD_1', 'WH_1',
'GDP_1', 'GDP_2', 'INT_1', 'IP_1',
'IM_1',
'TAX',
'M2CDN',
'CG', 'IG', 'IH', 'ST', 'EX',
'GDPDEF'
)]
Exog_df <- subset(Exog_df, Exog_df$YEAR >= 1982)
上記で抜き出した外生変数のテーブルには、外生というよりも「先決内生変数」と言われるものが混じっていて、たとえば「1期前のGDP」であれば、前のサイクルで推定された内生変数であるGDPの値を、当期の外生変数のようにして使うというものです。
先決内生変数については、見間違えないように、必要なもの以外はデータをNAに置き換えておくことにします。
1期のラグ変数なら1982年だけ残し、2期のラグ変数なら1983年まで残しておく必要があります。
remove_future_data <- function(df, valname, lag){
colnum <- which(colnames(df)==valname)
df[(lag+1):nrow(df),colnum] <- rep(NA, (nrow(df)-lag))
return(df)
}
Exog_df <- remove_future_data(df=Exog_df, valname='IP_1', lag=1)
Exog_df <- remove_future_data(df=Exog_df, valname='IM_1', lag=1)
Exog_df <- remove_future_data(df=Exog_df, valname='YD_1', lag=1)
Exog_df <- remove_future_data(df=Exog_df, valname='INT_1', lag=1)
Exog_df <- remove_future_data(df=Exog_df, valname='GDP_1', lag=1)
Exog_df <- remove_future_data(df=Exog_df, valname='GDP_2', lag=2)
内生変数の値を得る
「外生変数を与えると、連立方程式を解いて内生変数を返してくれる」という処理を書いていきます。
連立方程式を解くnleqslv()の使い方は次の通り。
未知数が並んだベクトルだけを引数に取る関数（f0とする）を定義して、nleqslv()に、関数foと「未知数の初期値ベクトル」を与えると、未知数の値をちょっとずつ変えてf0に与えながら、f0がゼロになる未知数の組み合わせを探してくれます。
つまりf0は、連立方程式が解けたらゼロベクトルになるように記述しておく必要があります。
たとえば、
INT = α + β1*ln(GDP) + β2*ln(M2CDN)
というような関係を想定しているのであれば、移項して、
INT -(α + β1*ln(GDP) + β2*ln(M2CDN))
というように記述しておくことになります。
で、以下のように定義する関数で、1期分の外生変数を与えると、その期の内生変数が得られます。
xの並び順に注意が必要です。
get_endo <- function(init, exo, year){
e <- subset(exo, exo$YEAR==year)
fn0 <- function(x){
output <- c(
x[1] - CP_fn(e$YD_1, e$WH_1),
x[2] - IP_fn(e$GDP_1, e$GDP_2, e$INT_1, e$IP_1),
x[3] - IM_fn(e$IM_1, x[6]),
x[4] - YD_fn(x[6], e$TAX),
x[5] - INT_fn(x[6], e$M2CDN),
x[6] - GDP_fn(x[1], x[2], e$CG, e$IG, e$IH, e$ST, e$EX, x[3]),
x[7] - GDPN_fn(x[6], e$GDPDEF)
)
return(output)
}
solution <- nleqslv(x=init, fn=fn0, method = "Newton")
return(solution$x)
}
内生変数を保存するテーブルを作っておきます。
Endog_df <- data.frame(YEAR = Exog_df$YEAR,
CP = rep(NA, nrow(Exog_df)),
IP = rep(NA, nrow(Exog_df)),
IM = rep(NA, nrow(Exog_df)),
YD = rep(NA, nrow(Exog_df)),
INT = rep(NA, nrow(Exog_df)),
GDP = rep(NA, nrow(Exog_df)),
GDPN = rep(NA, nrow(Exog_df))
)
次に、先ほど定義した関数を繰り返し使って各期の内生変数を順次求めていき、外生変数テーブルと内生変数テーブルをともに1期ずつ更新していく関数を書いておきます。
simulate <- function(init, exo_df, endo_df){
for (i in 1:nrow(exo_df)){
exo <- exo_df[i,]
sol <- get_endo(init=init, exo=exo_df, year=exo_df$YEAR[i])
endo_df[i,'CP'] <- sol[1]
endo_df[i,'IP'] <- sol[2]
endo_df[i,'IM'] <- sol[3]
endo_df[i,'YD'] <- sol[4]
endo_df[i,'INT'] <- sol[5]
endo_df[i,'GDP'] <- sol[6]
endo_df[i,'GDPN'] <- sol[7]
if(i<nrow(exo_df)){
exo_df[i+1,'IP_1'] <- sol[2]
exo_df[i+1,'IM_1'] <- sol[3]
exo_df[i+1,'YD_1'] <- sol[4]
exo_df[i+1,'INT_1'] <- sol[5]
exo_df[i+1,'GDP_1'] <- sol[6]
if(i>=2){
exo_df[i+1,'GDP_2'] <- endo_df[(i-1),'GDP']
}
}
}
output <- list(endo_df, exo_df)
return(output)
}
ではいよいよ下記でシミュレーションを実行します。
内生変数の初期値は、全部0とかにしててもきちんと推計できましたが、推計に使った元データの平均値を取って入れるようにしておきます。
initial <- c(mean(d_macro$CP), mean(d_macro$IP), mean(d_macro$IM), mean(d_macro$YD),
mean(d_macro$INT), mean(d_macro$GDP), mean(d_macro$GDPN))
results <- simulate(init=initial, exo_df=Exog_df, endo_df=Endog_df)
Endo_results <- results[[1]]
Exo_results <- results[[2]]
write.csv(Endo_results, 'results/Endo_results.csv', quote=FALSE, row.names=FALSE)
write.csv(Exo_results, 'results/Exo_results.csv', quote=FALSE, row.names=FALSE)
得られた結果を推計に使った元データと比較してみます。
GDPとINT（利子率）の予測値と実績値を並べておきます。実線が実績値、破線が予測値（シミュレーション値）です。
xlim=c(1982,2009)
ylim=c(300000, 600000)
plot(x=1982:2009, d_macro$GDP, type='l', lty=1,
xlim=xlim, ylim=ylim, ylab="GDP")
par(new=TRUE)
plot(x=1982:2009, Endo_results$GDP, type='l', lty=2,
xlim=xlim, ylim=ylim, axes = FALSE, xlab="", ylab="")
xlim=c(1982,2009)
ylim=c(0, 10)
plot(x=1982:2009, d_macro$INT, type='l', lty=1,
xlim=xlim, ylim=ylim, ylab="INT")
par(new=TRUE)
plot(x=1982:2009, Endo_results$INT, type='l', lty=2,
xlim=xlim, ylim=ylim, axes = FALSE, xlab="", ylab="")
教科書に従ってOLSで推定しましたが、2SLSにしても大して手間はかからないと思います。
