Verify SGU sampling of the last 4 parameters

[jlperla]

There may not be any issues at all, but the density plots of the last 4 parameters of SGU (1st order) are especially far off - especially for dynare. We want to make sure nothing loaded in the wrong order/etc.

We should verify that there aren't any major issues here. Some of the things to verify are

1. The ugliest ones are the rho, rho_u, rho_v, psi. As we know, the rho are the hardest to esimate so it might not even have any issues, just want to ensure that is the case.
2. The priors are correctly centered/etc. and identical between the models. In particular, we want to make sure we don't have an accidental trasncription error on the order of some parameter somewhere.
3. The pseudotrue is correctly defined for the DGP for both the dynare and the julia code
4. The sampling in dynare code is storing off the values correctly (i.e. https://github.com/HighDimensionalEconLab/HMCExamples.jl/blob/main/scripts/run_dynare_samplers/sgu_1.m#L10-L13) which is storing things in the order we think it is, and consistently.
   • Keep in mind the x2 is stored without any headers. I don't think this is causing trouble, but you never know where a reordering might be occuring in dyanre? THere is a variable called parameter_names in one of the dynare output files which seems to be in this order?
5. We convert this into a chain with https://github.com/HighDimensionalEconLab/HMCExamples.jl/blob/main/scripts/run_dynare_samplers/convert_dynare_output.jl#L5 etc. This takes in a list of parameters in order - which seems to match https://github.com/HighDimensionalEconLab/HMCExamples.jl/blob/main/scripts/run_dynare_samplers/convert_dynare_output.jl#L69

As I said, I can't find anything wrong which would lead to things loading the chain in the wrong order. But we should be dilligent before we declare dynare especeially horrible. If anything is wrong, my prior is on the priors and either there definition or connection to a DGP.

[donskerclass]

First investigating the most direct possibility, I can verify the ordering of the parameter names is correct by rerunning Dynare and checking that the results are indeed saved to CSV in the right order, and the parameter names in the Julia are in the order saved.

Priors do differ in the way reported in the paper: in Turing we use truncated normal priors for gamma and psi, while in Dynare we don't truncate. I forget exactly why we decided we were okay with this. I know Dynare only has limited options, but I think they are replicable in Turing if we need them. Maybe this was just something about avoiding negative values where the model would become undefined? Anyway, I doubt this is the source of large discrepancy, since little prior or posterior mass is in the truncated regions, but it's one of several differences.

The prior for beta_draw also differs due to rounding error in the gamma reparameterization: to be identical to Dynare it should be 4/9, 9, but 4/9 gets rounded to 0.444. I seriously doubt this makes a noticeable difference. I also double checked the reparameterization formulas and they look fine, at least if Wikipedia and the Distributions.jl docs are both correct.

The pseudotrues for generating the date are identical*, but note that we only have one sample of the data and we are generating different realizations for each, so unless the time length is long, we should expect estimates to differ just due to sampling variability. This is my number one most plausible explanation for posterior differences being not just possible but expected, in a way which does not show up for RBC where we use the same data. We could convert to using the same data with a few one-period timing shifts, and if this were the start of the project as opposed to the end I probably would have said we should do that, for comparability.

*We have a tiny bit of rounding error in the observation noise since Dynare uses variance and Julia uses standard deviation, but it's about 0.15 percent relative difference so likely negligible.

Also in the different generations Dynare starts from the ergodic instead of the steady state: this again would not matter with a very long data set, but will in finite samples. This probably makes a small but systematic difference: my guess is it's swamped by the sampling variability, but would show up at least a bit if we did frequentist diagnostics for this model, god forbid. Note that the samplers are informed of this difference, which is a major reason we were okay with different data, since it is not easy/possible to get the other samplers working with the different initial condition, so neither is misspecified for the own-model generated data, though the specifications do differ.

So, several minor differences are there, but no obvious definition or translation issues. The numbers are even fairly close for the first order. The second order a bit less so, but given sample quality that's plausible.

My main view on this is there's no reason, given the obvious source of variability that we're using two completely different series realizations to do the estimation, to think the differences are primarily due to errors.

There are things we could do to mitigate differences, if we were so inclined:

The big project one would be: remove model differences like the prior truncation, rounding, and especially the initial condition, simulate an identical data source, shift the columns for the Dynare to match the timing convention, and then estimate both on the same data. This would then allow more direct comparisons. This is way too much work and I don't recommend it.

A minor but feasible change would be to just plot the SGU estimates on different graphs, to make clear that comparison is not direct. This is problematic because we're short for space, but an option.

Or we can just leave as is and flag all the differences in the paper and footnotes, which we do now. That's my preferred option.

[jlperla]

just leave as is and flag all the differences in the paper and footnotes, which we do now. That's my preferred option.

@donskerclass Ditto. Correct is essential, but I see no reason to try to align these things any further if we can explain the differences. The models have subtle differences so it might not even work out. If referees are confused then we can put on separate graphs.

To me the reason to keep on the same graph is it gives a sense of the horizonal scale for how wonky the densities look like - which is our main purpose here. We can just add a comment on why the posterior modes are likely in different locations and it still servces that purpose.

So if there is nothing else to do, I think you can close this issue?