TCN DSGE example done

Examples/DSGE/SNM-TCN/DSGE_Example.jl → Examples/DSGE/SNM-TCN/Estimate.jl

@@ -1,5 +1,5 @@
 ## This does TCN neural net estimation for the DSGE example
-using PrettyTables, Pkg, DelimitedFiles, Distributions, LinearAlgebra
+using PrettyTables, Pkg, DelimitedFiles, Distributions, LinearAlgebra, MCMCChains, StatsPlots
 cd(@__DIR__)
 Pkg.activate(".")
 # defines the net and the DSGE model, and needed functions
@@ -39,7 +39,7 @@ X[:, 1, :] = Float32.(data)
 
 # compute mean and cov of moments, for obj fn and proposal
 function simmomentscov(θ::Vector{Float64}, S::Int64)
-m = Float64.(UntransformParameters(net(MakeData(θ, covreps, CKmodel)))')
+m = Float64.(UntransformParameters(net(MakeData(θ, S, CKmodel)))')
 mean(m, dims=1)[:], cov(m)
 end
 
@@ -107,10 +107,10 @@ end
 # proposal
 covreps = 1000
 _,Σₚ = simmomentscov(θnn, covreps)
-δ = 1.0 # tuning
+δ = 0.75 # tuning
 
 # define objective and proposal
-S = 20  # number of simulations for moments
+S = 50  # number of simulations for moments
 obj = θ -> bmsmobjective(θ, θnn, S)
 prop = θ -> proposal(θ, δ, Σₚ)
 
@@ -121,7 +121,7 @@ chn = Chains(chain[:,1:end-2], ["β", "γ", "ρ₁", "σ₁", "ρ₂", "σ₂",
 plot(chn)
 savefig("chain.png")
 display(chn)
-pretty_table([θtrue θnn mean(chain[:,1:end-2],dims=1)[:]], header = (["θtrue", "θnn", "θmcmc"]))
+pretty_table([TrueParameters() θnn mean(chain[:,1:end-2],dims=1)[:]], header = (["θtrue", "θnn", "θmcmc"]))
 
 end
 main()

Examples/DSGE/SNM-TCN/EstimateOnGPU.jl

@@ -0,0 +1,123 @@
+# try doing MCMC on GPU
+# this appears to be about 40% faster than on CPU, rough guess
+
+using PrettyTables, Pkg, DelimitedFiles, Distributions, LinearAlgebra, MCMCChains, StatsPlots
+cd(@__DIR__)
+Pkg.activate(".")
+# defines the net and the DSGE model, and needed functions
+include("Setup.jl")
+
+function main()
+
+net = load_trained()
+Flux.testmode!(net)
+net |> gpu
+CKmodel |> gpu
+
+## Now, let's move on to Bayesian MSM using either the typical data set, or generate a new one
+# load the data
+data = readdlm("dsgedata.txt")
+# transform the data the same way as was used to train net
+data .-= [0.84, 0.69, 0.33, 0.05, 1.72]'
+data ./= [0.51, 0.44, 0.36, 0.018, 0.34]'
+X = zeros(Float32, 160, 1, 5)
+X[:, 1, :] = Float32.(data)
+data = tabular2conv(permutedims(Float32.(X), (3, 2, 1)))
+data |> gpu
+fit = net(data)
+## This is the raw TCN estimate using the official data set
+θnn = Float64.(UntransformParameters(fit))[:]
+θnn |> gpu
+
+
+#################### Define functions for MCMC ###############################
+
+# compute mean and cov of moments, for obj fn and proposal
+function simmomentscov(θ::Vector{Float64}, S::Int64)
+    data = MakeData(θ, S, CKmodel)
+    fit = net(data)
+    m = Float64.(UntransformParameters(fit)')
+    mean(m, dims=1)[:], cov(m)
+end
+
+# CUE objective, written to MAXIMIZE
+@inbounds function bmsmobjective(θ::Vector{Float64}, θnn::Vector{Float64}, S::Int64)    # Make sure the trial parameter value is in the support
+    InSupport(θ) || return -Inf
+    # Compute simulated moments and covariance
+    θbar, Σ = simmomentscov(θ, S)
+    n = 160 # sample size
+    Σ *= n * (1+1/S) # 1 for θhat, 1/S for θbar
+    isposdef(Σ) || return -Inf
+    err = sqrt(n)*(θnn-θbar) 
+    W = inv(Σ)
+    -0.5*dot(err, W, err)
+end
+
+# proposal: MVN random walk
+@inbounds function proposal(current::Vector{Float64}, δ::Float64, Σ::Array{Float64})
+    rand(MvNormal(current, δ*Σ))
+end
+
+@views function mcmc(
+    θ::Vector{Float64}; # TODO: prior? not needed at present, as priors are uniform
+    Lₙ::Function, proposal::Function, burnin::Int=100, N::Int=1_000,
+    verbosity::Int=10
+)
+    Lₙθ = Lₙ(θ) # Objective at data moments value
+    naccept = 0 # Number of acceptance / rejections
+    accept = false
+    acceptance_rate = 1f0
+    chain = zeros(N, size(θ, 1) + 2)
+    for i ∈ 1:burnin+N
+        θᵗ = proposal(θ) # new trial value
+        Lₙθᵗ = Lₙ(θᵗ) # Objective at trial value
+        # Accept / reject trial value
+        accept = rand() < exp(Lₙθᵗ - Lₙθ)
+        if accept
+            # Replace values
+            θ = θᵗ
+            Lₙθ = Lₙθᵗ
+            # Increment number of accepted values
+            naccept += 1
+        end
+        # Add to chain if burnin is passed
+        # @info "current log-L" Lₙθ
+        if i > burnin
+            chain[i-burnin,:] = vcat(θ, accept, Lₙθ)
+        end
+        # Report
+        if verbosity > 0 && mod(i, verbosity) == 0
+            acceptance_rate = naccept / verbosity
+            @info "Current parameters (iteration i=$i)" round.(θ, digits=3)' acceptance_rate
+            naccept = 0
+        end
+    end
+    return chain
+end
+
+#################### End Define functions for MCMC ###############################
+
+
+## set up proposal and chain
+
+# proposal
+covreps = 1000
+_,Σₚ = simmomentscov(θnn, covreps)
+δ = 0.75 # tuning
+
+# define objective and proposal
+S = 50  # number of simulations for moments
+obj = θ -> bmsmobjective(θ, θnn, S)
+prop = θ -> proposal(θ, δ, Σₚ)
+
+## run the chain
+chain = mcmc(θnn, Lₙ=obj, proposal=prop, N=2000)
+# report results
+chn = Chains(chain[:,1:end-2], ["β", "γ", "ρ₁", "σ₁", "ρ₂", "σ₂", "nss"])
+plot(chn)
+savefig("chain.png")
+display(chn)
+pretty_table([TrueParameters() θnn mean(chain[:,1:end-2],dims=1)[:]], header = (["θtrue", "θnn", "θmcmc"]))
+
+end
+main()

Examples/DSGE/SimulatedNeuralMoments/Estimate.jl

@@ -141,6 +141,6 @@ chain = SimulatedNeuralMoments.mcmc(start, finallength, lnL, model, nnmodel, nni
 chn = Chains(chain[:,1:end-2], ["β", "γ", "ρ₁", "σ₁", "ρ₂", "σ₂", "nss"])
 plot(chn)
 #savefig("chain.png") # saved one used 5000 draws, current settings are fewer
-pretty_table([θtrue θnn mean(chain[:,1:end-2],dims=1)[:]], header = (["θtrue", "θnn", "θmcmc"]))
+pretty_table([TrueParameters() θnn mean(chain[:,1:end-2],dims=1)[:]], header = (["θtrue", "θnn", "θmcmc"]))
 display(chn)
 

econometrics.lyx

@@ -53873,6 +53873,32 @@ the 95% confidence intervals always contain the true parameter values
 \end_inset
 
 
+\end_layout
+
+\begin_layout Standard
+The chain follows.
+ It seems to mix reasonably well,
+ but in practice,
+ one would like a longer chain to reduce the variance of the posterior,
+ especially in the tails,
+ in order to get reliable confidence intervals.
+\end_layout
+
+\begin_layout Standard
+\begin_inset Graphics
+	filename Examples/DSGE/SimulatedNeuralMoments/chain.png
+	width 15cm
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset Newpage newpage
+\end_inset
+
+
 \end_layout
 
 \begin_layout Subsection
@@ -53921,30 +53947,48 @@ This approach has been found to have lower RMSE than the maximum likelihood esti
 \begin_layout Itemize
 This paper is certainly not the last word on the topic,
  but it does show that neural nets have much promise for improving the efficiency of simulated method of moments type estimators.
+\end_layout
+
+\end_deeper
+\begin_layout Itemize
+software to train a net,
+ and to do estimation is at 
+\begin_inset ERT
+status open
+
+\begin_layout Plain Layout
+
+
+\backslash
+href{./Examples/DSGE/SNM-TCN/Estimate.jl}{Estimate.jl}
+\end_layout
+
+\end_inset
+
+ and associated files.
 \begin_inset Newpage newpage
 \end_inset
 
 
 \end_layout
 
-\end_deeper
 \begin_layout Standard
 Monte Carlo results for 1000 replications of the TCN neural net estimator,
- for the simple DSGE model are
+ for the simple DSGE model are below.
 \end_layout
 
 \begin_layout Standard
 \begin_inset Graphics
 	filename Examples/DSGE/SNM-TCN/montecarlo.png
-	width 20cm
+	width 15cm
 
 \end_inset
 
 
 \end_layout
 
-\begin_layout Itemize
-if you compare these with the SNM neural moments that are based on a summary statistic (
+\begin_layout Standard
+Below are the results using the SNM neural net that takes a summary statistic as the input (
 \begin_inset CommandInset ref
 LatexCommand ref
 reference "SNM for DSGE results"
@@ -53956,7 +54000,25 @@ nolink "false"
 \end_inset
 
 ),
- you will see that these are considerably better for 
+ for comparison.
+\begin_inset Newline newline
+\end_inset
+
+
+\begin_inset Graphics
+	filename Examples/DSGE/SimulatedNeuralMoments/nnfit.png
+	width 15cm
+	special height=15cm
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Itemize
+the TCN results are better overall,
+ both in terms of bias and rmse,
+ and are considerably better for 
 \begin_inset Formula $\gamma$
 \end_inset
 
@@ -53965,16 +54027,62 @@ nolink "false"
 
 \begin_layout Itemize
 This shows that there may be information lost when using summary statistics,
- and that we can do better with a well-chosen neural net.
+ and that we can do better with a well-chosen neural net that uses the full sample information.
+ However,
+ more evidence would be needed to confirm these results.
 \end_layout
 
 \begin_layout Itemize
+95% confidence intervals are tight,
+ and contain the true parameter values,
+ for all parameters.
+\begin_inset Newpage newpage
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+Results for CUE-GMM using the TCN fit as the simulated moments are
+\end_layout
+
+\begin_layout Standard
+\begin_inset Graphics
+	filename Examples/DSGE/SNM-TCN/Results.png
+	width 20cm
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Itemize
+The point estimates are very close to the true values,
+ and the standard deviations are small,
+ and are,
+ in general,
+ better than what we have seen with other estimators.
+
+\begin_inset Newpage newpage
+\end_inset
+
+
+\end_layout
+
+\begin_layout Itemize
+The chains follow.
+ They look pretty good.
+ It would be better to use longer chains,
+ to smooth out the bumps in the marginal densities.
+\end_layout
+
+\begin_layout Standard
+\begin_inset Graphics
+	filename Examples/DSGE/SNM-TCN/chain.png
+	width 20cm
+
+\end_inset
 
-\size footnotesize
-Note:
- Still need to add Bayesian MSM,
- using the neural net fit as the moments for MSM.
-\size default
 
 \begin_inset Newpage newpage
 \end_inset