Commit 2369cc4e authored by Thomas Brand's avatar Thomas Brand

Initial commit.

parents
---
output:
html_document:
code_folding: hide
toc: true
bibliography: fs2000.bib
---
```{r, message=FALSE, warning=FALSE, echo=FALSE, results='hide'}
rm(list = ls())
if (!"pacman" %in% installed.packages()[,"Package"]) install.packages("pacman", repos='http://cran.r-project.org')
pacman::p_load(curl,dplyr,magrittr,tidyr,ggplot2,lubridate,knitr,rsdmx,zoo,R.matlab,stringr)
opts_chunk$set(message=FALSE, warning=FALSE, cache=FALSE, comment=NA)
# theme for ggplot
theme <- theme_bw()+ theme(strip.background=element_blank(),
strip.text=element_text(size=15),
title=element_text(size=16),
panel.border=element_blank(),
panel.grid.major=element_line(size=1),
legend.key=element_rect(colour="white"),
legend.position="bottom",
legend.text=element_text(size=10),
axis.text=element_text(size=10))
blueObsMacro <- "#0D5BA4"
matlab_path <- '/usr/local/bin/matlab'
dynare_path <- '~/Documents/dynare/matlab'
#dynare_path <- '/usr/lib/dynare/matlab'
```
This post replicates the estimation of the cash in advance model (termed M1 in the paper) described in @Scho00, with updated data. This implementation was written by Michel Juillard and is given in the `examples` folder of Dynare.
# Data
The data are retrieved directly through the Widukind databases and are stocked in `fs2000_data.csv`. They are therefore different from the original dataset used by @Scho00.
```{r}
url_widukind <- "http://widukind-api.cepremap.org/api/v1/sdmx"
url_provider <- "/BEA/data"
url_dataset_series <- "/nipa-section1-10106-q/A191RX1.Q"
# Real Gross Domestic Product, unit : billions 2009 $, SA annual rate
url <- paste0(url_widukind,url_provider,url_dataset_series)
data_sdmx <- readSDMX(url)
df_gdp <- data.frame(data_sdmx, var="gdp")
url_dataset_series <- "/nipa-section1-10109-q/A191RD3.Q"
# Gross Domestic Product: Implicit Price Deflator, SA
url <- paste0(url_widukind,url_provider,url_dataset_series)
data_sdmx <- readSDMX(url)
df_defgdp <- data.frame(data_sdmx, var="defgdp")
bea_data <-
rbind(df_gdp,df_defgdp) %>%
select(var, TIME_PERIOD, OBS_VALUE) %>%
rename(time = TIME_PERIOD,
values = OBS_VALUE) %>%
mutate(time=as.Date(as.yearqtr(time,format="%YQ%q")),
values=as.numeric(values))
url_provider <- "/OECD/data"
url_dataset_series <- "/MEI/USA.LFWA64TT.STSA.Q"
# unit: 1000 person (source: OECD, MEI)
url <- paste0(url_widukind,url_provider,url_dataset_series)
data_sdmx <- readSDMX(url)
df_pop <- data.frame(data_sdmx, var="pop")
pop <-
df_pop %>%
select(var, TIME_PERIOD, OBS_VALUE) %>%
rename(time = TIME_PERIOD,
values = OBS_VALUE) %>%
mutate(time=as.Date(as.yearqtr(time,format="%Y-Q%q")),
values=as.numeric(values))
rawdata <- rbind(bea_data,pop)
fs2000_data <-
rawdata %>%
spread(key = var, value = values) %>%
transmute(time=time,
gdp_rpc=1e+9*gdp/pop,
defgdp = defgdp,
gy_obs = gdp_rpc/lag(gdp_rpc),
gp_obs = defgdp/lag(defgdp))
fs2000_data %>%
transmute(time=gsub(" ","",as.yearqtr(time)),
gy_obs,
gp_obs) %>%
na.omit() %>%
write.csv("fs2000_data.csv", row.names=FALSE)
```
# Model
The equations are taken from @Naso94. Note that there is an initial minus sign missing in equation (A1), p. S63. The prior distribution follows the one originally specified in Schorfheide's paper. Note that the $\beta$ prior for $\rho$ implies an asymptote and corresponding prior mode for $\rho$ at 0. It is generally recommended to avoid this extreme type of prior.
## Preamble
```{r}
cat(readLines("fs2000.mod")[1:12], sep= '\n')
```
## Model equations
```{r}
cat(readLines("fs2000.mod")[14:29], sep= '\n')
```
## Steady state
```{r}
cat(readLines("fs2000.mod")[31:58], sep= '\n')
```
## Shocks
```{r}
cat(readLines("fs2000.mod")[60:63], sep= '\n')
```
## Estimation
```{r}
cat(readLines("fs2000.mod")[65:77], sep= '\n')
```
## Computation
```{r}
cat(readLines("fs2000.mod")[79], sep= '\n')
```
You can run Dynare directly from this file, but before that you must check the `dynare_path` and `matlab_path` and you must have a version of Dynare >= 4.4 (in order to read csv file for estimation).
```{r eval=F}
system(paste(matlab_path, '-nosplash -nodisplay -r "addpath', dynare_path , '; dynare fs2000; exit" '))
```
# Results
```{r}
listPriorShape <- list("Beta"=1,
"Gamma"=2,
"Gaussian"=3,
"Inverse Gamma 1"=4,
"Uniform"=5,
"Inverse Gamma 2"=6,
"Weibull"=8)
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
```{r}
# prior
dname <- unlist(subset(M_,rownames(M_)=="dname"))
param_estim <- as.vector(unlist(mode[["parameter_names"]]))
priorshape <- readMat(paste0("./",dname,"/prior/definition.mat"))[[1]][[1]][,1]
priormean <- unlist(subset(oo_,rownames(oo_)=="prior")[[1]][1])
priorvar <- diag(matrix(unlist(subset(oo_,rownames(oo_)=="prior")[[1]][3]),nrow=length(priormean),byrow=TRUE))
dataPrior <- data.frame(names = param_estim,
prior_shape = priorshape,
prior_mean = priormean,
prior_stdv = sqrt(priorvar))
dataPrior$prior_shape <- factor(dataPrior$prior_shape)
levels(dataPrior$prior_shape) <- listPriorShape
# post
postmode <- subset(oo_,row.names(oo_)=="posterior_mode")
postnames <- c(row.names(postmode[[1]][[1]]),
row.names(postmode[[1]][[2]]))
dataPost <- data.frame(names=postnames,
unlist(subset(oo_,row.names(oo_)=="posterior_mode")[[1]]),
unlist(subset(oo_,row.names(oo_)=="posterior_mean")[[1]]),
unlist(subset(oo_,row.names(oo_)=="posterior_hpdinf")[[1]]),
unlist(subset(oo_,row.names(oo_)=="posterior_hpdsup")[[1]]),
unlist(subset(oo_,row.names(oo_)=="posterior_std")[[1]]))
priorPost<-merge(dataPrior,dataPost,by="names")
colnames(priorPost)<-c('parameter','Prior density','Prior mean','Prior stdv',
"Post. mode","Post. mean",
"Prob. Interval 10%","Prob. Interval 90%",
"Post. stdv")
kable(priorPost)
```
## Bayesian impulse response function
```{r}
```
# References
\ No newline at end of file
% This file was created with JabRef 2.10.
% Encoding: ISO8859_1
@Article{Naso94,
Title = {Testing the implications of long-run neutrality for monetary business cycle models},
Author = {Nason, James M and Cogley, Timothy},
Journal = {Journal of Applied Econometrics},
Year = {1994},
Number = {S1},
Pages = {S37--S70},
Volume = {9},
Doi = {10.1002/jae.3950090504}
}
@Article{Scho00,
Title = {{Loss function-based evaluation of DSGE models}},
Author = {Schorfheide, Frank},
Journal = {Journal of Applied Econometrics},
Year = {2000},
Number = {6},
Pages = {645--670},
Volume = {15},
Doi = {10.1002/jae.582},
File = {Scho00.pdf:Scho00.pdf:PDF},
ISSN = {1099-1255},
Owner = {tbrand},
Timestamp = {2016.06.29}
}
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m;
parameters alp bet gam mst rho psi del;
alp = 0.33;
bet = 0.99;
gam = 0.003;
mst = 1.011;
rho = 0.7;
psi = 0.787;
del = 0.02;
model;
dA = exp(gam+e_a);
log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
-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);
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;
end;
steady_state_model;
dA = exp(gam);
gst = 1/dA;
m = mst;
khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
n = xist/(nust+xist);
P = xist + nust;
k = khst*n;
l = psi*mst*n/( (1-psi)*(1-n) );
c = mst/P;
d = l - mst + 1;
y = k^alp*n^(1-alp)*gst^alp;
R = mst/bet;
W = l/n;
ist = y-c;
q = 1 - d;
e = 1;
gp_obs = m/dA;
gy_obs = dA;
end;
steady;
check;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
estimated_params;
alp, beta_pdf, 0.356, 0.02;
bet, beta_pdf, 0.993, 0.002;
gam, normal_pdf, 0.0085, 0.003;
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;
end;
varobs gp_obs gy_obs;
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;
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