Commit d7497786 authored by Thomas Brand's avatar Thomas Brand

Add bayesian irf analysis.

parent 2369cc4e
......@@ -83,13 +83,13 @@ fs2000_data <-
transmute(time=time,
gdp_rpc=1e+9*gdp/pop,
defgdp = defgdp,
gy_obs = gdp_rpc/lag(gdp_rpc),
gp_obs = defgdp/lag(defgdp))
gyObs = gdp_rpc/lag(gdp_rpc),
gpObs = defgdp/lag(defgdp))
fs2000_data %>%
transmute(time=gsub(" ","",as.yearqtr(time)),
gy_obs,
gp_obs) %>%
gyObs,
gpObs) %>%
na.omit() %>%
write.csv("fs2000_data.csv", row.names=FALSE)
```
......@@ -149,47 +149,14 @@ listPriorShape <- list("Beta"=1,
results <- readMat("fs2000_results.mat",fixNames = F)
mode <- readMat("fs2000_mode.mat",fixNames = F)
options_ <- results[["options_"]]
M_ <- results[["M_"]]
oo_ <- results[["oo_"]]
dataset_ <- results[["dataset_"]]
varendo_code <- unlist(subset(M_,row.names(M_)=="endo_names")) %>%
str_trim()
varendo_name <- unlist(subset(M_,row.names(M_)=="endo_names_long")) %>%
str_trim()
varendo <- data.frame(code=varendo_code,
name=varendo_name,
stringsAsFactors = F)
varexo_code <- unlist(subset(M_,row.names(M_)=="exo_names")) %>%
str_trim()
varexo_name <- unlist(subset(M_,row.names(M_)=="exo_names_long")) %>%
str_trim()
varexo <- data.frame(code=varexo_code,
name=varexo_name,
stringsAsFactors=F)
param_code <- unlist(subset(M_,row.names(M_)=="param_names")) %>%
str_trim()
param_name <- unlist(subset(M_,row.names(M_)=="param_names_long")) %>%
str_trim()
param <- data.frame(code=param_code,
name=param_name,
stringsAsFactors=F)
var_code_name <- rbind(varendo,varexo,param)
varobs <- unlist(subset(options_,row.names(options_)=="varobs"))
time <- as.data.frame(dataset_$date$time) %>%
unite(time,c(V1,V2),sep="Q")
time <- as.character(as.Date(as.yearqtr(time$time)))
```
## Posterior for the parameters of the structural model
Our results look very similar to the posterior means and standard errors for the parameters of the structural model (M1) of @Scho00 (Table II, p.659).
```{r}
# prior
dname <- unlist(subset(M_,rownames(M_)=="dname"))
......@@ -224,16 +191,63 @@ colnames(priorPost)<-c('parameter','Prior density','Prior mean','Prior stdv',
"Post. mode","Post. mean",
"Prob. Interval 10%","Prob. Interval 90%",
"Post. stdv")
priorPost[,-c(1,2)]<-round(priorPost[,-c(1,2)],4)
kable(priorPost)
```
## Bayesian impulse response function
The shape of our impulse response function are similar to the posterior distribution of impulse response functions of @Scho00 (Figure 1, first row, p.664).
```{r}
birf <- subset(oo_,row.names(oo_)=="PosteriorIRF")[[1]][[1]]
names <- as.vector(row.names(subset(birf,row.names(birf)=="Mean")[[1]]))
meanbirf <- data.frame(meas="mean",
var = names,
matrix(unlist(subset(birf,row.names(birf)=="Mean")),
nrow=length(names),
byrow=T))
infbirf <- data.frame(meas="inf",
var = names,
matrix(unlist(subset(birf,row.names(birf)=="HPDinf")),
nrow=length(names),
byrow=T))
supbirf <- data.frame(meas="sup",
var = names,
matrix(unlist(subset(birf,row.names(birf)=="HPDsup")),
nrow=length(names),
byrow=T))
birf <- rbind(meanbirf,infbirf,supbirf) %>%
separate(col=var,into=c("var","shock"))
colnames(birf) <- gsub("X","",colnames(birf))
birf_df <-
birf %>%
gather(time,values,-meas,-shock,-var) %>%
mutate(time=as.numeric(time),
values=100*values) %>%
spread(meas,values)
listVar <- list("1- Output growth" = "gyObs",
"2- Inflation" = "gpObs")
birf_df$var <- factor(birf_df$var)
levels(birf_df$var)<-listVar
```
```{r fig.align="center", fig.height=4}
birf_df %>% filter(shock=="eM") %>%
ggplot(aes(x=time,y=mean,group=shock))+
geom_line(size=0.8,colour="#1F77B4")+
geom_ribbon(aes(ymin=inf,ymax=sup),alpha=0.6,fill="#D6E3F3")+
facet_wrap(~var,scales="free_y")+
xlab(NULL) + ylab("Percent")+theme_bw()+theme+
ggtitle("Monetary shock")
```
# References
\ No newline at end of file
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m;
var m P c e W R k d n l gyObs gpObs y dA;
varexo eA eM;
parameters alp bet gam mst rho psi del;
......@@ -12,20 +12,20 @@ psi = 0.787;
del = 0.02;
model;
dA = exp(gam+e_a);
log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
dA = exp(gam+eA);
log(m) = (1-rho)*log(mst) + rho*log(m(-1))+eM;
-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
W = l/n;
-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
R = P*(1-alp)*exp(-alp*(gam+eA))*k(-1)^alp*n^(-alp)/W;
1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+eA))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
c+k = exp(-alp*(gam+eA))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+eA))*k(-1);
P*c = m;
m-1+d = l;
e = exp(e_a);
y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
gy_obs = dA*y/y(-1);
gp_obs = (P/P(-1))*m(-1)/dA;
e = exp(eA);
y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+eA));
gyObs = dA*y/y(-1);
gpObs = (P/P(-1))*m(-1)/dA;
end;
steady_state_model;
......@@ -50,16 +50,16 @@ steady_state_model;
e = 1;
gp_obs = m/dA;
gy_obs = dA;
gpObs = m/dA;
gyObs = dA;
end;
steady;
check;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
var eA; stderr 0.014;
var eM; stderr 0.005;
end;
estimated_params;
......@@ -70,11 +70,11 @@ mst, normal_pdf, 1.0002, 0.007;
rho, beta_pdf, 0.129, 0.223;
psi, beta_pdf, 0.65, 0.05;
del, beta_pdf, 0.01, 0.005;
stderr e_a, inv_gamma_pdf, 0.035449, inf;
stderr e_m, inv_gamma_pdf, 0.008862, inf;
stderr eA, inv_gamma_pdf, 0.035449, inf;
stderr eM, inv_gamma_pdf, 0.008862, inf;
end;
varobs gp_obs gy_obs;
varobs gpObs gyObs;
estimation(order=1, datafile=fs2000_data, loglinear, mh_replic=2000, mh_nblocks=2, mh_jscale=0.8, bayesian_irf, irf=40, nograph) gy_obs gp_obs;
estimation(order=1, datafile=fs2000_data, loglinear, mh_replic=2000, mh_nblocks=2, mh_jscale=0.8, bayesian_irf, irf=40, nograph) gyObs gpObs;
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment