Skip to content

Commit

Permalink
Finished wrteup
Browse files Browse the repository at this point in the history
  • Loading branch information
@jlperla
jlperla committed Jun 2, 2023
1 parent 94f4d23 commit 36c2d7d
Show file tree
Hide file tree
Showing 3 changed files with 36 additions and 138 deletions.
39 changes: 13 additions & 26 deletions README.md
View file
Original file line number Diff line number Diff line change
@@ -1,13 +1,19 @@
# DSSM HMC Examples

This package provides the compete replication for [__"Differentiable State-Space Models and Hamiltonian Monte Carlo Estimation"__](https://www.jesseperla.com/publication/diff-state-space/diff-state-space.pdf) by David Childers, Jesus Fernandez-Villaverde, Jesse Perla, Christopher Rackauckas, and Peifan Wu.

The algorithms are mostly implemented in:
1. [DifferentiableStateSpaceModels.jl](https://github.com/HighDimensionalEconLab/DifferentiableStateSpaceModels.jl) for the differentiable perturbation method solver
2. [DifferenceEquations.jl](https://github.com/SciML/DifferenceEquations.jl/) for the differentiable solver for the Kalman filter, simulations of state space models, and likelihoods.

## Installation
Do `julia --project -e "using Pkg; Pkg.instantiate()"`
With Julia 1.9 and this repository cloned to your local machine, execute `julia --project -e "using Pkg; Pkg.instantiate()"`

## CLI Usage
To use with default options, and specifying a directory
```bash
julia --project --threads auto bin/fit_rbc_1_kalman.jl --results_path ./.results/rbc_1_kalman --overwrite_results true --num_samples 100
julia --project --threads auto bin/fit_rbc_1_joint.jl --results_path ./.results/rbc_1_joint --overwrite_results true --num_samples 100
julia --project bin/fit_rbc_1_kalman.jl --results_path ./.results/rbc_1_kalman --overwrite_results true --num_samples 100
julia --project bin/fit_rbc_1_joint.jl --results_path ./.results/rbc_1_joint --overwrite_results true --num_samples 100
```
(although the `--threads auto` may or may not be useful in this example.)

Expand All @@ -20,8 +26,6 @@ julia --project --threads auto bin/fit_rbc_1_kalman.jl --results_path ./.results
A few features for sampling and output:
- `--target_acceptance_rate 0.65` set the target acceptance for NUTS
- `--max_depth 5` set the max_depth for NUTS
- `--discard_initial 100` etc. discards draws as a warmup
- `--save_jls true` save the entire chain as an (unportable) JLS file. Default is true
- `--init_params_file data/my_file.csv` uses that file as the initial condition for sampling

In all cases, for long-burnins you can find the final draw of the chain as the `last_draw.csv` file, which can be used with the `init_params_file` argument after moving/renaming
Expand All @@ -32,31 +36,14 @@ julia --threads auto -e 'using Pkg; Pkg.add(\"PackageCompiler\")'
julia --project --threads auto ./deps/create_sysimage.jl
```

Grab coffee. This will take at least 15-30 minutes to run. It will help startup latency (maybe 2-3 minutes less to start sampling), but not as much as you would hope. Zygote.jl is not amenable to caching compilation.
Grab coffee. This will take at least 15-30 minutes to run. It will help startup latency (maybe 2-3 minutes less to start sampling), but not as much as you would hope. Zygote.jl is not amenable to caching compilation when used with Turing.

After: when you use vscode it will also load this custom sysimage as long as you have the `Julia: Use custom sysimage` option enabled.

On the commandline you will need to manually provide it manually. For example,
On the commandline you will need to provide the image manually. For example,
```bash
julia --project --sysimage JuliaSysimage.dll --threads auto bin/fit_rbc_1_kalman.jl --results_path ./.results/rbc_1_kalman --overwrite_results true --num_samples 1000
```

# Analyzing the Chain
## Loading files

To load a file into a chain with some path which saved with baseline serialization
```julia
using Serialization, HMCExamples, DelimitedFiles, MCMCChains
chain = open(deserialize, joinpath(pkgdir(HMCExamples), ".results/rbc_1_joint/chain.jls"))
```

## Accessing Chains
To extract the entire chain for some parameters
```
vals = get(chain, [:α,:β_draw])
```

Or to get a slice of all parameters for the last draw in the chain
```julia
last_draw = chain.value[end,:,1][chain.name_map.parameters] |> Array
```
## Paper Replication
To replicate all of the results in the paper, see the [README.md](scripts/README.md) in the scripts directory.
115 changes: 13 additions & 102 deletions scripts/README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ To run the primary experiments, execute:
bash scripts/run_samplers/baseline_experiments.sh
```

The other scripts are left separate to make running them in parallel easier. These run thousands of examples and may take **3-4 days to run** in parallel, depending on your computer. However, you can edit the shell scripts to run these in parallel with little effort.
The other scripts are left separate to make running them in parallel easier. These run thousands of examples and may take **3-4 days to run** in parallel, depending on your computer. However, you can execute in parallel on a multi-core machine to speed up the process.

```bash
bash scripts/run_samplers/rbc_1_joint_frequentist.sh
Expand Down Expand Up @@ -84,108 +84,19 @@ Then, you can run the following (in separate terminals if you wish, as the resul
cd scripts/run_dynare_samplers
matlab -nosplash -nodesktop -r "run('rbc_1_robustness.m');exit;"
matlab -nosplash -nodesktop -r "run('rbc_2_robustness.m');exit;"
cd ../.. # go back to the main directory
```

Finally, with all of the dynare results complete run
## Generating figures and tables
Assuming that you have either executed the above steps, or downloaded a preexisting `.replication_results` and put it local to your computer, you can generate all of the figures and tables to `.paper_results` by running `bash scripts/generate_paper_results.sh`. These will be placed in the `.paper_results` directory.

```bash
julia --project scripts/run_dynare_samplers/convert_dynare_output.jl
```

## Generating figures and tables
Assuming that you have either executed the above steps, or downloaded a `.replication_results` and put it local to your computer you can generate all of the figures and tables to `.paper_results` by runing `bash scripts/generate_paper_results.sh`. This will take a few hours to run. The individual scripts are:

1. `/generate_paper_results/baseline_tables.py` creates the primary tables for all of the experiments. Note that the `convert_dynare_output.jl` is necessary to execute beforehand.
2.


## TODO BELOW
## Producing summary tables and plots

On Ubuntu, Python or Python3 should be installed by default. If not, do `sudo apt update` followed by `sudo apt install python3` and `sudo apt install python3-pip`. If the command is present as `python3`, replace all instances of `python` in the text below with `python3`.\
Following this, run `cd HMCExamples.jl` followed by `pip install -r scripts/figplots/requirements.txt`.

Next, you will need to instantiate the Julia project associated with the scripts, which is separate from the main project file. Run `julia --project=scripts -e "using Pkg; Pkg.instantiate()"`.

In addition, you will need to have all experiments copied with the proper filepaths into a folder called `.experiments` in the `HMCExamples.jl` directory. If adding experiments manually, check the corresponding plot/table script for the filepaths being used.

Run `mkdir .tables .results .figures` to prevent filepath errors when running the scripts.

The following options are available:

### Dynare RBC
- Dynare Log Parser
- Use this to read the Dynare logs if you saved it to a file earlier.
- `python scripts/figplots/parse_dynare_log.py`

- Summary Statistics
- First, pull the Matlab chain files into Julia by running `julia --project=scripts scripts/figplots/dynare_to_julia.jl`.
- Next, parse the results into tables using `python scripts/figplots/sumstats_dynare.py`.

- Robustness Plots
- `julia --project=scripts --threads auto scripts/figplots/robustness_dynare.jl`
- Ideally, this command should be run on a machine with as many threads as possible due to the slow execution time of StatsPlots.

### Dynare SGU
- Summary Statistics
- First, pull the Matlab chain files into Julia by running `julia --project=scripts scripts/sgu_replication/dynare_to_julia.jl`.
- Next, parse the results into tables using `python scripts/sgu_replication/sumstats_dynare.py`.
- Plots are available under the Julia SGU section.

### Julia RBC
- Numerical Error Presence Validation
- `python scripts/figplots/scan_for_errors.py`

- Robustness Rhat Analysis/Validation
- `python scripts/figplots/scan_rhat.py`

- Summary Statistics
- `python scripts/figplots/sumstats_julia.py`

- Summary Plots
- `julia --project=scripts scripts/figplots/statplots_julia.jl`

- Platform Cross-Comparison, Scatterplot, Inferred Shocks
- `julia --project=scripts scripts/figplots/otherplots.jl`

- Robustness Plots
- `julia --project=scripts scripts/figplots/robustness_julia.jl`

- Frequentist Statistics
- First, run `julia --project=scripts scripts/figplots/frequentist_julia.jl`.
- After, run `python scripts/figplots/frequentist_julia.py`.

### Julia FVGQ
- Summary Statistics
- `python scripts/fvgq_replication/sumstats.py`

- Summary Plots
- `julia --project=scripts scripts/fvgq_replication/statplots.jl`

### Julia SGU
- Summary Statistics
- `python scripts/sgu_replication/sumstats_julia.py`

- Summary Plots
- `julia --project=scripts scripts/sgu_replication/statplots.jl`

### Julia RBC with Student's T Shocks
- Summary Statistics
- `python scripts/rbc_student_t_replication/sumstats_julia.py`

- Summary Plots
- `julia --project=scripts scripts/rbc_student_t_replication/statplots.jl`

### Julia RBC with Stochastic Volatility, 1st Order
- Summary Statistics
- `python scripts/rbc_volatility_replication/sumstats_julia.py`

- Summary Plots
- `julia --project=scripts scripts/rbc_volatility_replication/statplots.jl`

### Julia RBC with Stochastic Volatility, 2nd Order, SGU Form
- Summary Statistics
- `python scripts/rbc_sv_replication/sumstats_julia.py`

- Summary Plots
- `julia --project=scripts scripts/rbc_sv_replication/statplots.jl`
The individual scripts this calls are:
1. `convert_frequentist_output.jl`: Converts the output from the frequentist experiments into a format that can be used by the plotting scripts.
2. `convert_dynare_output.jl`: Converts the output from the dynare experiments into a format that can be used by the plotting scripts consistent with the Julia chains.
3. `baseline_figures.jl`: Generates all figures except for the RBC robustness examples
4. `rbc_robustness_figures.jl`: Generates the RBC robustness figures
5. `baseline_tables.py`: Generates all tables except for the RBC frequentist tables
6. `rbc_frequentist_tables.py`: Generates the RBC frequentist tables

Constants such as the number of samples are extracted from metadata in the `.replication_results` directory.
20 changes: 10 additions & 10 deletions scripts/generate_paper_results/baseline_figures.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,7 @@ function dynare_rbc_comparison(julia_small_run, julia_big_run, dynare_run, pseud
title!(density_plot_alpha, title[1])
title!(density_plot_beta, title[2])
title!(density_plot_rho, title[3])
# ylabel!(density_plot_alpha, "Sample Value")
# ylabel!(density_plot_alpha, "Density")
ylabel!(density_plot_beta, "")
ylabel!(density_plot_rho, "")

Expand Down Expand Up @@ -203,7 +203,7 @@ savefig(plt, ".paper_results/rbc_2_joint_200_density_traceplots.png")
# 2nd RBC dynare comparison
show_pseudo_true = false
density_plot_alpha, density_plot_beta, density_plot_rho = dynare_rbc_comparison("rbc_2_joint_200", "rbc_2_joint_200_long", "rbc_2_200_dynare", rbc_pseudotrue, title;show_pseudo_true, left_margin = 7mm, top_margin = 5mm)
plt = plot(density_plot_alpha, density_plot_beta, density_plot_rho; layout=(1,3), size = (1200, 300), legend=:bottomright, lw=2, ylabel=["Sample Value" "" ""])
plt = plot(density_plot_alpha, density_plot_beta, density_plot_rho; layout=(1,3), size = (1200, 300), legend=:bottomright, lw=2, ylabel=["Density" "" ""])
savefig(plt, ".paper_results/rbc_2_dynare_comparison.png")

# Scatterplots
Expand Down Expand Up @@ -255,17 +255,17 @@ sgu_vars_titles = [L"\alpha", L"\beta_{draw}", L"\gamma", L"\rho", L"\rho_u", L"
sgu_pseudotrues = [0.32 4 2.0 0.42 0.2 0.4 0.000742]

plt = dynare_sgu_comparison_1("sgu_1_kalman_200", "sgu_1_joint_200", "sgu_1_200_dynare", sgu_include_vars, sgu_pseudotrues, sgu_vars_titles;show_pseudo_true, layout=(4,2), size = (750, 800), legend=:topleft, left_margin = 1mm, rightmargin=3mm)
ylabel!(plt[1], "Sample Value")
ylabel!(plt[3], "Sample Value")
ylabel!(plt[5], "Sample Value")
ylabel!(plt[7], "Sample Value")
ylabel!(plt[1], "Density")
ylabel!(plt[3], "Density")
ylabel!(plt[5], "Density")
ylabel!(plt[7], "Density")
savefig(plt, ".paper_results/sgu_1_comparison.png")

plt = dynare_sgu_comparison_2("sgu_2_joint_200", "sgu_2_200_dynare", sgu_include_vars, sgu_pseudotrues, sgu_vars_titles;show_pseudo_true, layout=(4,2), size = (750, 800), legend=:topleft, left_margin = 1mm, rightmargin=3mm)
ylabel!(plt[1], "Sample Value")
ylabel!(plt[3], "Sample Value")
ylabel!(plt[5], "Sample Value")
ylabel!(plt[7], "Sample Value")
ylabel!(plt[1], "Density")
ylabel!(plt[3], "Density")
ylabel!(plt[5], "Density")
ylabel!(plt[7], "Density")
savefig(plt, ".paper_results/sgu_2_comparison.png")

# # Traceplots combined. Left out of paper for now
Expand Down

0 comments on commit 36c2d7d

Please sign in to comment.