now using Dynare for ML and MCMC of DSGE model

Examples/DSGE/Bayesian/CKmcmc.mod → Examples/DSGE/Bayesian/CKmcmc_cn.mod

@@ -1,21 +1,21 @@
-// The two shock model
+// The two shock model from Creel and Kristensen "Indirect Likelihood Inference"
 // This takes steady state of hours as given, and computes the psi
 // parameter that corresponds.
-close all;
 
 var y c k n invest z1 z2 MUC MUL r w;
 varexo e1 e2;
 
 parameters alppha betta delta gam  nss rho1 sigma1 rho2 sigma2 psi c1 iss yss kss css;
-set_param_value('alppha',0.33)
-set_param_value('betta', 0.99)
-set_param_value('delta',0.025)
-set_param_value('gam', 2.0)
-set_param_value('rho1', 0.9)
-set_param_value('sigma1', 0.02)
-set_param_value('rho2', 0.7)
-set_param_value('sigma2', 0.01)
-set_param_value('nss', 1/3)
+
+alppha = 0.33;
+betta = 0.99;
+delta = 0.025;
+gam = 2.0;
+rho1 = 0.9;
+sigma1 = 0.02;
+rho2 = 0.7;
+sigma2 = 0.01;
+nss =1/3;
 
 model;
   MUC = (c)^(-gam);
@@ -56,20 +56,27 @@ z1 = 0;
 z2 = 0;
 end;
 
+steady;
+
+stoch_simul(order=1);
+
 
 estimated_params;
- betta, uniform_pdf, , ,  0.95, 0.9999; 
- gam, uniform_pdf, , , 0.0,  5.0;
- rho1, uniform_pdf, , , 0.0,  1.0; 
- sigma1, uniform_pdf, , , 0.0,  0.1;
- rho2, uniform_pdf, , ,   0.0,  1.0; 
- sigma2, uniform_pdf, , , 0.0,  0.1;
- nss, uniform_pdf, , , 6/24, 9/24; 
+// start values for the optimization 
+// the numbers after the comment are the true 
+// values used to generate the data.
+betta   ,   0.99    ,   uniform_pdf     ,   0.9724  ,   0.01299;    // 0.99
+gam     ,   2.0     ,   uniform_pdf     ,   2.5     ,   1.4434;       // 2.0
+rho1    ,   0.9     ,   uniform_pdf     ,   0.4975  ,   0.28723;      // 0.9
+sigma1  ,   0.02    ,   uniform_pdf     ,   0.05    ,   0.02887;   // 0.02
+rho2    ,   0.7     ,   uniform_pdf     ,   0.4975  ,   0.28723;       // 0.7
+sigma2  ,   0.01    ,   uniform_pdf     ,   0.05    ,   0.02887;   // 0.01
+nss     ,   .3333   ,   uniform_pdf     ,   0.3125  ,    0.03608;      // 1/3
 end;
 
-varobs c n;  // experiment choosing one or two from y c n MPK MPL (c and n work well)
+varobs c n;  // experiment choosing one or two from y c n r w
 
-// short MCMC for illustration only 
-estimation(order=1, datafile=dsgedata, nobs=160, mh_replic=10000, mh_nblocks=1, mh_jscale=0.8, mh_init_scale=5, irf=40);
+// Does a single fairly short chain 
+estimation(order=1, datafile='dsgedata.csv', mh_replic=20000, mh_jscale=2.) ;
 
 

Examples/DSGE/Bayesian/CKmcmc_yw.mod

@@ -0,0 +1,82 @@
+// The two shock model from Creel and Kristensen "Indirect Likelihood Inference"
+// This takes steady state of hours as given, and computes the psi
+// parameter that corresponds.
+
+var y c k n invest z1 z2 MUC MUL r w;
+varexo e1 e2;
+
+parameters alppha betta delta gam  nss rho1 sigma1 rho2 sigma2 psi c1 iss yss kss css;
+
+alppha = 0.33;
+betta = 0.99;
+delta = 0.025;
+gam = 2.0;
+rho1 = 0.9;
+sigma1 = 0.02;
+rho2 = 0.7;
+sigma2 = 0.01;
+nss =1/3;
+
+model;
+  MUC = (c)^(-gam);
+  MUL = psi*exp(z2);
+  r = alppha * exp(z1) * (k(-1))^(alppha-1) * n^(1-alppha);
+  w = (1-alppha)*exp(z1)* (k(-1))^alppha * n^(-alppha);
+  MUC = betta*MUC(+1) * (1 + r(+1) - delta);
+  MUL/MUC = w;
+  z1 = rho1*z1(-1) + sigma1*e1;
+  z2 = rho2*z2(-1) + sigma2*e2;
+  y = exp(z1) * ((k(-1))^alppha) * (n^(1-alppha));
+  invest = y - c;
+  k = invest + (1-delta)*k(-1);
+end;
+
+shocks;
+var e1 = 1;
+var e2 = 1;
+end;
+
+steady_state_model;
+c1 = ((1/betta + delta - 1)/alppha)^(1/(1-alppha));
+kss = nss/c1;
+iss = delta*kss;
+yss = kss^alppha * nss^(1-alppha);
+css = yss - iss;
+psi =  (css^(-gam)) * (1-alppha) * (kss^alppha) * (nss^(-alppha));
+k = kss;
+n = nss;
+c = css;
+y = yss;
+invest = iss;
+MUC = c^(-gam);
+r = alppha  * k^(alppha-1) * n^(1-alppha);
+w = (1-alppha)* (k)^alppha * n^(-alppha);
+MUL = w*MUC;
+z1 = 0;
+z2 = 0;
+end;
+
+steady;
+
+stoch_simul(order=1);
+
+
+estimated_params;
+// start values for the optimization 
+// the numbers after the comment are the true 
+// values used to generate the data.
+betta   ,   0.99    ,   uniform_pdf     ,   0.9724  ,   0.01299;    // 0.99
+gam     ,   2.0     ,   uniform_pdf     ,   2.5     ,   1.4434;       // 2.0
+rho1    ,   0.9     ,   uniform_pdf     ,   0.4975  ,   0.28723;      // 0.9
+sigma1  ,   0.02    ,   uniform_pdf     ,   0.05    ,   0.02887;   // 0.02
+rho2    ,   0.7     ,   uniform_pdf     ,   0.4975  ,   0.28723;       // 0.7
+sigma2  ,   0.01    ,   uniform_pdf     ,   0.05    ,   0.02887;   // 0.01
+nss     ,   .3333   ,   uniform_pdf     ,   0.3125  ,    0.03608;      // 1/3
+end;
+
+varobs y w;  // experiment choosing one or two from y c n r w
+
+// Does a single fairly short chain 
+estimation(order=1, datafile='dsgedata.csv', mh_replic=20000, mh_jscale=7.) ;
+
+

Examples/DSGE/Bayesian/cleanup

@@ -1,7 +1,5 @@
 #!/bin/bash
-mv CKmcmc.mod temp
-rm -R -f +CKmcmc
-rm -R -f CK*
-rm *.mat
-rm *.dat
-mv temp CKmcmc.mod
+rm -R -f CKmcmc_cn
+rm -R -f CKmcmc_yw
+rm *.log
+

Examples/DSGE/Bayesian/dsgedata.csv

@@ -0,0 +1 @@
+../GenData/dsgedata.csv

Examples/DSGE/GenData/CKNlib.jl

@@ -0,0 +1,152 @@
+using Econometrics, Statistics, LinearAlgebra, SolveDSGE
+
+# this block reads and processes the file, leave it be
+process_model("CK.txt")
+const dsge = retrieve_processed_model("CK_processed.txt")
+
+# solve model and simulate data
+function dgp(θ, dsge, reps, rndseed=1234)
+    model_params, ss, noise = ParamsAndSS(θ)
+    dsge = assign_parameters(dsge, model_params)
+    scheme = PerturbationScheme(ss, 1.0, "third")
+    solution = solve_model(dsge, scheme)
+    burnin = 500
+    nobs = 160
+    data = simulate(solution, ss[1:3], reps*(burnin+nobs); seed = rndseed)
+    # add measurement error, where st. dev. is a given proportion of steady state
+    data[4:8,:] += randn(5, size(data,2)) .* noise .* ss[4:8]
+    # the next returns reps data sets, in an array of arrays
+    data = [data[4:8, (nobs+burnin)*i-(nobs+burnin)+1+burnin:i*(nobs+burnin)]' for i = 1:reps]
+end
+
+function bad_data(data)
+    any(isnan.(data)) || any(isinf.(data)) || any(std(data,dims=1) .==0.0) || any(data .< 0.0)
+end
+
+# this gives a vector of vectors, each a statistic drawn at the parameter value
+function auxstat(θ, reps)
+    auxstat.(dgp(θ, dsge, reps, rand(1:Int64(1e12))))
+end
+
+# These are the candidate auxiliary statistics for ABC estimation of
+# the simple DSGE model of Creel and Kristensen (2013)
+@views function auxstat(data)
+    # check for nan, inf, no variation, or negative, all are reasons to reject
+    if bad_data(data)
+        return zeros(39)
+    else    
+        # known parameters
+        α = 0.33
+        δ = 0.025
+        # recover capital (see notes)
+        hours = data[:,3]
+        intrate = data[:,4]
+        wages = data[:,5]
+        capital = α /(1.0-α )*hours.*wages./intrate
+        # treat all variables
+        logdata = log.([data capital])
+        # logs
+        logoutput = logdata[:,1];   # output
+        logcons = logdata[:,2];     # consumption
+        loghours = logdata[:,3];    # hours
+        logintrate = logdata[:,4];  # intrate
+        logwages = logdata[:,5];    # wages
+        logcapital = logdata[:,6]
+        # rho1, sig1
+        e = logoutput-α*logcapital-(1.0-α)*loghours 
+        y = e[2:end]
+        x = e[1:end-1]
+        rho1 = x\y
+        u = y-x*rho1
+        sig1 = sqrt(mean(u .^2))
+        # gam, rho2, sig2 (1/MRS=wage)
+        x = [ones(size(logcons,1)) logcons]
+        b = x\logwages
+        e = logwages-x*b
+        y = e[2:end]
+        x = e[1:end-1]
+        rho2 = x\y
+        u = y-x*rho2
+        sig2 = sqrt(mean(u .^2))
+        # standard devs. and correlations
+        m = mean(logdata, dims=1)
+        d = (logdata .- m)  # keep means and std. devs., the VAR uses standardized and normalized   
+        # AR(1)
+        maxlag = 1
+        y = d[2:end,:]
+        x = d[1:end-1,:]
+        n = size(y,1)
+        rhos = zeros(6)
+        es = zeros(n,6)
+        for i = 1:6
+            rho = x[:,i]\y[:,i]
+            rhos[i] = rho
+            es[:,i] = y[:,i]-rho.*x[:,i]
+        end        
+        Z = vcat(rho1, sig1, b, rho2, sig2, m[:], rhos, vech(cov(es)))
+    end
+    Z
+end
+
+function TrueParameters()
+ [
+ 0.99,  # β
+ 2.0,   # γ     
+ 0.9,   # ρ₁  
+ 0.02,  # σ₁   
+ 0.7,   # ρ₂  
+ 0.01,  # σ₂   
+ 8.0/24.0,  # nss
+ 0.02, # measurement error has std. dev. 2% of steady state
+ 0.02,
+ 0.02,
+ 0.02,
+ 0.02]
+end    
+
+# measurement error std. dev. is between 0 and 10 percent of steady state
+function PriorSupport()
+    lb = [0.95, 0.0, 0.0, 0.0, 0.0, 0.0, 6.0/24.0, 0.,0.,0.,0.,0.]
+    ub = [0.995, 5.0, 0.995, 0.1, 0.995, 0.1, 9.0/24.0,0.1,0.1,0.1,0.1,0.1]
+    lb,ub
+end    
+
+function PriorDraw()
+    lb, ub = PriorSupport()
+    θ = (ub-lb).*rand(size(lb,1)) + lb
+end    
+
+function InSupport(θ)
+    lb,ub = PriorSupport()
+    all(θ .>= lb) & all(θ .<= ub)
+end
+
+# in auxstats, we already code bad data, not needed here
+function GoodData(z)
+    true
+end
+
+function Prior(θ)
+    InSupport(θ) ? 1.0 : 0.0
+end    
+
+function ParamsAndSS(θ)
+    α = 0.33
+    δ = 0.025
+    β, γ, ρ₁, σ₁, ρ₂, σ₂, nss, σy, σc, σn, σw, σr = θ
+    c1 = ((1/β  + δ - 1)/α)^(1/(1-α))
+    kss = nss/c1
+    iss = δ*kss
+    yss = kss^α * nss^(1-α)
+    css = yss - iss;
+    MUCss = css^(-γ)
+    rss = α * kss^(α-1) * nss^(1-α)
+    wss = (1-α)* (kss)^α * nss^(-α)
+    MULss = wss*MUCss
+    ψ =  (css^(-γ)) * (1-α) * (kss^α) * (nss^(-α))
+    model_params = [β, γ, ρ₁ , σ₁, ρ₂, σ₂, ψ]
+    noise_params =  [σy, σc, σn, σw, σr] 
+    ss = [0.0, 0.0, kss, yss, css, nss, rss, wss, MUCss, MULss]
+    return model_params, ss, noise_params
+end   
+

Examples/DSGE/GenData/GenNoisyData.jl

@@ -0,0 +1,24 @@
+using SolveDSGE, Plots, DelimitedFiles, DataFrames, CSV
+cd(@__DIR__)
+
+# use the library with measurement error
+include("CKNlib.jl")
+
+if !isfile("./CK_processed.txt") 
+    process_model("./CK.txt")
+end    
+global const dsge = retrieve_processed_model("./CK_processed.txt")
+
+function GenData(silent=false)
+# this block reads and processes the file, leave it be
+data = dgp(TrueParameters(), dsge, 1)[1]
+df = DataFrame(data, ["output", "cons", "hours","r","w"])
+if !silent
+    display(plot(data, legend=:outertopright, label=["output" "cons" "hours" "r" "w"]))
+    savefig("dsgedata.svg")
+    display(df)
+end    
+#writedlm("dsgedata.txt", data)
+return true
+end
+