text:
"""
```
function load_posterior_moments(m; load_bands = true, include_fixed = false)
```

Load posterior moments (mean, std) of parameters for a particular sample, and optionally
also load 5% and 95% lower and upper bands.

### Keyword Arguments
- `cloud::ParticleCloud`: Optionally pass in a cloud that you want to load the sample from. If the cloud is non-empty then the model object will only be used to find fixed indices and parameter tex labels
- `load_bands::Bool`: Optionally include the 5% and 95% percentiles for the sample of parameters in the returned df
- `include_fixed::Bool`: Optionally include the fixed parameters in the returned df
- `excl_list::Vector{Symbol}`: List parameters by their key that you want to exclude from
loading

### Outputs
- `df`: A dataframe containing the aforementioned moments/bands
"""
function load_posterior_moments(m::AbstractDSGEModel;
                                cloud::Union{SMC.Cloud,DSGE.Cloud,ParticleCloud} = ParticleCloud(m, 0),
                                load_bands::Bool = true,
                                include_fixed::Bool = false,
                                excl_list::Vector{Symbol} = Vector{Symbol}(undef, 0),
                                weighted::Bool = true)
    parameters = m.parameters

    # Read in Posterior Draws
    if cloud_isempty(cloud)
        if get_setting(m, :sampling_method) == :MH
            params = load_draws(m, :full)
            params = get_setting(m, :sampling_method) == :MH ? thin_mh_draws(m, params) : params # TODO
            params = params'
            weights = Vector{Float64}(undef, 0)
            weighted = false
        elseif get_setting(m, :sampling_method) == :SMC
            cloud = get_cloud(m)
            params = get_vals(cloud)
            weights = get_weights(cloud)
        else
            @error "Invalid sampling method"
        end
    else
        params = get_vals(cloud)
        weights = get_weights(cloud)
   end

    # Index out the fixed parameters
    if include_fixed
        free_indices   = 1:n_parameters(m)
    else
        free_indices   = findall(x -> x.fixed == false, parameters)
        parameters = parameters[free_indices]
        params = params[free_indices, :]
    end

    # Remove excluded parameters
    included_indices = setdiff(1:length(free_indices), calculate_excluded_indices(m, excl_list, include_fixed = include_fixed))
    parameters = parameters[included_indices]
    params = params[included_indices, :]
    tex_labels = [DSGE.detexify(parameters[i].tex_label) for i in 1:length(parameters)]

    load_posterior_moments(params, weights, tex_labels, weighted = weighted, load_bands = load_bands)
end

# Aggregating many ParticleClouds to calculate moments across
# multiple SMC estimations
function load_posterior_moments(m::AbstractDSGEModel,
                                clouds::Union{Vector{ParticleCloud},Vector{SMC.Cloud}};
                                load_bands::Bool = true,
                                include_fixed::Bool = false,
                                excl_list::Vector{Symbol} = Vector{Symbol}(undef, 0),
                                weighted::Bool = true)
    n_params = n_parameters(m)
    n_particles = length(clouds[1])
    n_clouds = length(clouds)
    n_free_params = length(Base.filter(x -> x.fixed == false, m.parameters))

    p_mean = Array{Float64}(undef, n_free_params, n_clouds)
    p_lb = Array{Float64}(undef, n_free_params, n_clouds)
    p_ub = Array{Float64}(undef, n_free_params, n_clouds)

    for (i, c) in enumerate(clouds)
        df_i = load_posterior_moments(m; weighted = weighted, cloud = c, load_bands = load_bands, include_fixed = include_fixed, excl_list = excl_list)
        p_mean[:, i] = df_i[!,:post_mean]
        p_lb[:, i] = df_i[!,:post_lb]
        p_ub[:, i] = df_i[!,:post_ub]
    end
    param = load_posterior_moments(m; cloud = clouds[1], load_bands = true, include_fixed = false)[!,:param]
    df = DataFrame(param = param, post_mean = dropdims(mean(p_mean, dims = 2), dims = 2), post_std = dropdims(std(p_mean, dims = 2),dims = 2), post_lb = dropdims(mean(p_lb, dims = 2), dims = 2), post_ub = dropdims(mean(p_ub, dims = 2), dims = 2))

    return df
end

# Base method
# Assumes n_parameters x n_draws
function load_posterior_moments(params::Matrix{Float64}, weights::Vector{Float64}, tex_labels::Vector{String}; weighted::Bool = true, load_bands::Bool = true)
    if size(params, 1) > size(params, 2)
        @warn "`params` argument to load_posterior_moments seems to be oriented incorrectly.
        The argument should be n_params x n_draws. Currently, size(params) = ($(size(params, 1)), $(size(params, 2)))"
    end
    if weighted
        params_mean = vec(mean(params, Weights(weights), dims = 2))
        params_std = vec(std(params, Weights(weights), 2, corrected = false))
    else
        params_mean = vec(mean(params, dims = 2))
        params_std  = vec(std(params, dims = 2))
    end

    df = DataFrame()
    df[!,:param] = tex_labels
    df[!,:post_mean] = params_mean
    df[!,:post_std]  = params_std

    if load_bands
        post_lb = Vector{Float64}(undef, length(params_mean))
        post_ub = similar(post_lb)
        for i in 1:length(params_mean)
            if weighted
                post_lb[i] = quantile(params[i, :], Weights(weights), .05)
                post_ub[i] = quantile(params[i, :], Weights(weights), .95)
            else
                post_lb[i] = quantile(params[i, :], .05)
                post_ub[i] = quantile(params[i, :], .95)
            end
        end
    end

    df[!,:post_lb] = post_lb
    df[!,:post_ub] = post_ub

    return df
end

# For calculating the indices to exclude when making posterior interval plots
function calculate_excluded_indices(m::AbstractDSGEModel, excl_list::Vector{Symbol}; include_fixed::Bool = false)
    parameters = deepcopy(m.parameters)

    # If fixed parameters are excluded, the indices need to be calculated
    # with respect to only the free parameters
    if include_fixed
        return findall(x -> x in excl_list, [parameters[i].key for i in 1:length(parameters)])
    else
        # Remove the fixed parameters
        Base.filter!(x -> x.fixed == false, parameters)
        return findall(x -> x in excl_list, [parameters[i].key for i in 1:length(parameters)])
    end
end

"""
```
moment_tables(m; percent = 0.90, subset_inds = 1:0, subset_string = "",
    groupings = Dict{String, Vector{Parameter}}(), use_mode = false,
    tables = [:prior_posterior_means, :moments, :prior, :posterior],
    caption = true, outdir = "", verbose = :none)
```

Computes prior and posterior parameter moments. Tabulates prior mean, posterior
mean, and bands in various LaTeX tables. These tables will be saved in `outdir`
if it is nonempty, or else in `tablespath(m, \"estimate\")`.

### Inputs

- `m::AbstractDSGEModel`: model object

### Keyword Arguments

- `percent::AbstractFloat`: the percentage of the mass of draws from
  Metropolis-Hastings included between the bands displayed in output tables.
- `subset_inds::AbstractRange{Int64}`: indices specifying the draws we want to use
- `subset_string::String`: short string identifying the subset to be
  appended to the output filenames. If `subset_inds` is nonempty but
  `subset_string` is empty, an error is thrown
- `groupings::Dict{String, Vector{Parameter}}`: see `?parameter_groupings`
- `use_mode::Bool`: use the modal parameters instead of the mean in the
  prior_posterior_means table
- `tables::Vector{Symbol}`: which tables to produce
- `caption::Bool`: whether to include table captions
- `outdir::String`: where to save output tables
- `verbose::Symbol`: desired frequency of function progress messages printed to
  standard out. One of `:none`, `:low`, or `:high`
"""
function moment_tables(m::AbstractDSGEModel; percent::AbstractFloat = 0.90,
                       subset_inds::AbstractRange{Int64} = 1:0, subset_string::String = "",
                       groupings::AbstractDict{String, Vector{Parameter}} = Dict{String, Vector{Parameter}}(),
                       tables = [:prior_posterior_means, :moments, :prior, :posterior],
                       caption = true, outdir = "",
                       verbose::Symbol = :low, use_mode::Bool = false)

    ### 1. Load parameter draws from Metropolis-Hastings

    params = if !isempty(subset_inds)
        # Use subset of draws
        if isempty(subset_string)
            error("Must supply a nonempty subset_string if subset_inds is nonempty")
        end
        load_draws(m, :subset; subset_inds = subset_inds, verbose = verbose)
    else
        # Use all draws
        load_draws(m, :full; verbose = verbose)
    end

    ### 2. Compute posterior moments

    if use_mode
        post_mode = h5read(get_forecast_input_file(m, :mode), "params")
    end
    post_means = vec(mean(params, dims = 1))

    # Save posterior means
    basename = "paramsmean"
    if !isempty(subset_string)
        basename *= "_sub=$(subset_string)"
    end
    filename = workpath(m, "estimate", "$basename.h5")
    h5open(filename, "w") do file
        write(file, "post_means", post_means)
    end
    post_bands = permutedims(find_density_bands(params, percent; minimize = true))

    ### 3. Produce TeX tables

    if :prior_posterior_means in tables
        prior_posterior_table(m, use_mode ? post_mode : post_means;
                              subset_string = subset_string, groupings = groupings,
                              use_mode = use_mode, caption = caption, outdir = outdir)
    end

    if :moments in tables
        prior_posterior_moments_table(m, post_means, post_bands; percent = percent,
                                      subset_string = subset_string, groupings = groupings,
                                      caption = caption, outdir = outdir)
    end

    if :prior in tables
        prior_table(m, groupings = groupings, caption = caption, outdir = outdir)
    end

    if :posterior in tables
        posterior_table(m, post_means, post_bands, percent = percent,
                        subset_string = subset_string, groupings = groupings,
                        caption = caption, outdir = outdir)
    end
    #=
    if :mean_mode_moments in tables
        mean_mode_moments_table(m, post_means, post_bands; percent = percent,
                                subset_string = subset_string, groupings = groupings,
                                caption = caption, outdir = outdir)
    end
    =#
    println(verbose, :low, "Tables are saved as " * tablespath(m, "estimate", "*.tex"))
end

"""
```
moments(θ::Parameter)
```

If θ's prior is a `RootInverseGamma`, τ and ν. Otherwise, returns the mean
and standard deviation of the prior. If θ is fixed, returns `(θ.value, 0.0)`.
"""
function moments(θ::Parameter)
    if θ.fixed
        return θ.value, 0.0
    else
        prior = get(θ.prior)
        if isa(prior, RootInverseGamma)
            return prior.τ, prior.ν
        else
            return mean(prior), std(prior)
        end
    end
end

"""
```
prior_table(m; subset_string = "", groupings = Dict{String, Vector{Parameter}}(),
    caption = true, outdir = "")
```
"""
function prior_table(m::AbstractDSGEModel; subset_string::String = "",
                     caption = true, outdir = "",
             groupings::AbstractDict{String, Vector{Parameter}} = Dict{String, Vector{Parameter}}())

    if isempty(groupings)
        sorted_parameters = sort(m.parameters, by = (x -> x.key))
        groupings[""] = sorted_parameters
    end

    # Open the TeX file
    basename = "priors"
    if !isempty(subset_string)
        basename *= "_sub=$(subset_string)"
    end
    outfile = tablespath(m, "estimate", "$basename.tex")
    if !isempty(outdir)
        outfile = replace(outfile, dirname(outfile) => outdir)
    end
    fid = open(outfile, "w")

    # Write header
    write_table_preamble(fid)

    @printf fid "\\renewcommand*\\footnoterule{}"
    @printf fid "\\vspace*{.5cm}\n"
    @printf fid "{\\small\n"
    @printf fid "\\begin{longtable}{rlrr@{\\hspace{1in}}rlrr}\n"
    if caption
        @printf fid "\\caption{Priors}\n"
    end
    @printf fid "\\label{tab:param-priors}\n"
    @printf fid "\\\\ \\hline\n"

    @printf fid "& Dist & Mean & Std Dev & & Dist & Mean & Std Dev \\\\ \\hline\n"
    @printf fid "\\endhead\n"

    @printf fid "\\hline \\\\\n"
    @printf fid "\\multicolumn{8}{c}{\\footnotesize Note: For Inverse Gamma prior mean and SD, \$\\tau\$ and \$\\nu\$ reported.}\n"
    @printf fid "\\endfoot\n"

    # Map prior distributions to identifying strings
    distid(::Distributions.Uniform) = "Uniform"
    distid(::Distributions.Beta)    = "Beta"
    distid(::Distributions.Gamma)   = "Gamma"
    distid(::Distributions.Normal)  = "Normal"
    distid(::RootInverseGamma)      = "InvG"

    # Write priors
    for group_desc in keys(groupings)
        params = groupings[group_desc]

        # Take out anticipated shock SDs 2 to k - these priors are all the same
        antshock_params = [m[k] for k in [Symbol("σ_r_m$i") for i = 2:n_anticipated_shocks(m)]]
        params = setdiff(params, antshock_params)

        n_params = length(params)
        n_rows = convert(Int, ceil(n_params/2))

        # Write grouping description if not empty
        if !isempty(group_desc)
            @printf fid "\\multicolumn{8}{l}{\\textit{%s}} \\\\[3pt]\n" group_desc
        end

        # Write footnote about standard deviations of anticipated policy shocks
        function anticipated_shock_footnote(θ::Parameter)
            if n_anticipated_shocks(m) > 0 && θ.key == :σ_r_m1
                nantpad          = n_anticipated_shocks_padding(m)
                all_sigmas       = [m[Symbol("σ_r_m$i")]::Parameter for i = 1:nantpad]
                nonzero_sigmas   = Base.filter(θ -> !(θ.fixed && θ.value == 0), all_sigmas)
                n_nonzero_sigmas = length(nonzero_sigmas)

                if n_nonzero_sigmas > 1
                    text = "\$\\sigma_{ant1}\$ through \$\\sigma_{ant$(n_nonzero_sigmas)}\$ all have the same distribution."
                    @printf fid "\\let\\thefootnote\\relax\\footnote{\\centering %s}" text
                end
            end
        end

        for i = 1:n_rows
            # Write left column
            θ = params[i]
            (prior_mean, prior_std) = moments(θ)
            @printf fid "\$%s\$ &" θ.tex_label
            @printf fid " %s &" (θ.fixed ? "-" : distid(get(θ.prior)))
            @printf fid " %0.2f &" prior_mean
            if θ.fixed
                @printf fid " \\scriptsize{fixed} &"
            else
                @printf fid " %0.2f &" prior_std
            end
            anticipated_shock_footnote(θ)

            # Write right column if it exists
            if n_rows + i <= n_params
                θ = params[n_rows + i]
                (prior_mean, prior_std) = moments(θ)
                @printf fid " \$%s\$ &" θ.tex_label
                @printf fid " %s &" (θ.fixed ? "-" : distid(get(θ.prior)))
                @printf fid " %0.2f &" prior_mean
                if θ.fixed
                    @printf fid " \\scriptsize{fixed}"
                else
                    @printf fid " %0.2f" prior_std
                end
                anticipated_shock_footnote(θ)
            else
                @printf fid  "& & &"
            end

            # Add padding after last row in a grouping
            if i == n_rows
                @printf fid " \\\\[3pt]\n"
            else
                @printf fid " \\\\\n"
            end
        end
    end

    # Write footer
    write_table_postamble(fid; small = true)

    # Close file
    close(fid)
end

"""
```
posterior_table(m, post_means, post_bands; percent = 0.9, subset_string = "",
    groupings = Dict{String, Vector{Parameter}}(), caption = true, outdir = "")
```
"""
function posterior_table(m::AbstractDSGEModel, post_means::Vector, post_bands::Matrix;
                         percent::AbstractFloat = 0.9,
                         subset_string::String = "",
                         groupings::AbstractDict{String, Vector{Parameter}} = Dict{String, Vector{Parameter}}(),
                         caption::Bool = true,
                         outdir::String = "")

    if isempty(groupings)
        sorted_parameters = sort(m.parameters, by = (x -> x.key))
        groupings[""] = sorted_parameters
    end

    # Open the TeX file
    basename = "posterior"
    if !isempty(subset_string)
        basename *= "_sub=$(subset_string)"
    end
    outfile = tablespath(m, "estimate", "$basename.tex")
    if !isempty(outdir)
        outfile = replace(outfile, dirname(outfile) => outdir)
    end
    fid = open(outfile, "w")

    # Write header
    write_table_preamble(fid)

    @printf fid "\\renewcommand*\\footnoterule{}"
    @printf fid "\\vspace*{.5cm}\n"
    @printf fid "{\\small\n"
    @printf fid "\\begin{longtable}{rrc@{\\hspace{1in}}rrc}\n"
    if caption
        @printf fid "\\caption{Posteriors}\n"
    end
    @printf fid "\\label{tab:param-posteriors}\n"
    @printf fid "\\\\ \\hline\n"

    lb = (1 - percent)/2 * 100
    ub = 100 - lb
    @printf fid "& Mean & (p%0.0f, p%0.0f) & & Mean & (p%0.0f, p%0.0f) \\\\ \\hline\n" lb ub lb ub
    @printf fid "\\endhead\n"

    @printf fid "\\hline \\\\\n"
    @printf fid "\\endfoot\n"

    # Write priors
    for group_desc in keys(groupings)
        params = groupings[group_desc]
        n_params = length(params)
        n_rows = convert(Int, ceil(n_params/2))

        # Write grouping description if not empty
        if !isempty(group_desc)
            @printf fid "\\multicolumn{6}{l}{\\textit{%s}} \\\\[3pt]\n" group_desc
        end

        for i = 1:n_rows
            # Write left column
            θ = params[i]
            j = m.keys[θ.key]
            @printf fid "\$%s\$ &" θ.tex_label
            @printf fid " %0.2f &" post_means[j]
            if θ.fixed
                @printf fid " \\scriptsize{fixed} &"
            else
                @printf fid " (%0.2f, %0.2f) &" post_bands[j, :]...
            end

            # Write right column if it exists
            if n_rows + i <= n_params
                θ = params[n_rows + i]
                j = m.keys[θ.key]
                (prior_mean, prior_std) = moments(θ)
                @printf fid " \$%s\$ &" θ.tex_label
                @printf fid " %0.2f &" post_means[j]
                if θ.fixed
                    @printf fid " \\scriptsize{fixed}"
                else
                    @printf fid " (%0.2f, %0.2f)" post_bands[j, :]...
                end
            else
                @printf fid " & &"
            end

            # Add padding after last row in a grouping
            if i == n_rows
                @printf fid " \\\\[3pt]\n"
            else
                @printf fid " \\\\\n"
            end
        end
    end

    # Write footer
    write_table_postamble(fid; small = true)

    # Close file
    close(fid)
end

"""
```
prior_posterior_moments_table(m, post_means, post_bands; percent = 0.9,
    subset_string = "", groupings = Dict{String, Vector{Parameter}}(),
    caption = true, outdir = "")
```

Produces a table of prior means, prior standard deviations, posterior means, and
90% bands for posterior draws.
"""
function prior_posterior_moments_table(m::AbstractDSGEModel,
                                       post_means::Vector, post_bands::Matrix;
                                       percent::AbstractFloat = 0.9,
                                       subset_string::String = "",
                                       caption = true, outdir = "",
                 groupings::AbstractDict{String, Vector{Parameter}} = Dict{String, Vector{Parameter}}())

    if isempty(groupings)
        sorted_parameters = sort(m.parameters, by = (x -> x.key))
        groupings[""] = sorted_parameters
    end

    # Open the TeX file
    basename = "moments"
    if !isempty(subset_string)
        basename *= "_sub=$(subset_string)"
    end
    outfile = tablespath(m, "estimate", "$basename.tex")
    if !isempty(outdir)
        outfile = replace(outfile, dirname(outfile) => outdir)
    end
    fid = open(outfile, "w")

    # Write header
    write_table_preamble(fid)

    @printf fid "\\vspace*{.5cm}\n"
    @printf fid "{\\small\n"
    @printf fid "\\begin{longtable}{lcccccc}\n"
    if caption
        @printf fid "\\caption{Parameter Estimates}\n"
    end
    @printf fid "\\\\ \\hline\n"

    # Two-row column names. First row is multicolumn, where entries (i,str) are `i` columns
    # with content `str`.
    colnames0 = [(1,""), (3,"Prior"),(3,"Posterior")]
    colnames = ["Parameter", "Type", "Mean", "SD", "Mean",
                "$(100*percent)\\% {\\tiny Lower Band}",
                "$(100*percent)\\% {\\tiny Upper Band}"]

    @printf fid "\\multicolumn{%d}{c}{%s}" colnames0[1][1] colnames0[1][2]
    for col in colnames0[2:end]
        @printf fid " & \\multicolumn{%d}{c}{%s}" col[1] col[2]
    end
    @printf fid " \\\\\n"
    @printf fid "%s" colnames[1]
    for col in colnames[2:end]
        @printf fid " & %s" col
    end
    @printf fid " \\\\\n"
    @printf fid "\\cmidrule(lr){1-1} \\cmidrule(lr){2-4} \\cmidrule(lr){5-7}\n"
    @printf fid "\\endhead\n"

    @printf fid "\\hline\n"
    @printf fid "\\\\ \\multicolumn{7}{c}{\\footnotesize Note: For Inverse Gamma (IG) prior mean and SD, \$\\tau\$ and \$\\nu\$ reported.}\n"
    @printf fid "\\endfoot\n"

    # Map prior distributions to identifying strings
    distid(::Distributions.Uniform) = "U"
    distid(::Distributions.Beta)    = "B"
    distid(::Distributions.Gamma)   = "G"
    distid(::Distributions.Normal)  = "N"
    distid(::RootInverseGamma) = "IG"

    # Write parameter moments
    # sorted_parameters = sort(m.parameters, by = (x -> x.key))
    for group_desc in keys(groupings)
        params = groupings[group_desc]

        # Write grouping description if not empty
        if !isempty(group_desc)
            @printf fid "\\multicolumn{7}{c}{\\textit{%s}} \\\\[3pt]\n" group_desc
        end

        for param in params
            index = m.keys[param.key]
            (prior_mean, prior_std) = moments(param)

            @printf fid "\$%4.99s\$ & " param.tex_label
            @printf fid "%s & " (param.fixed ? "-" : distid(get(param.prior)))
            @printf fid "%8.3f & " prior_mean
            @printf fid "%8.3f & " prior_std
            @printf fid "%8.3f & " post_means[index]
            @printf fid "%8.3f & %8.3f \\\\\n" post_bands[index, :]...
        end
    end

    # Write footer
    write_table_postamble(fid; small=true)

    # Close file
    close(fid)
end

"""
```
prior_posterior_table(m, post_values; subset_string = "",
    groupings = Dict{String, Vector{Parameter}}(), use_mode = false,
    caption = true, outdir = "")
```

Produce a table of prior means and posterior means or mode.
"""
function prior_posterior_table(m::AbstractDSGEModel, post_values::Vector;
                               subset_string::String = "",
                               groupings::AbstractDict{String,Vector{Parameter}} = Dict{String, Vector{Parameter}}(),
                               caption = true, outdir = "",
                               use_mode::Bool = false)

    if isempty(groupings)
        sorted_parameters = sort(m.parameters, by = (x -> x.key))
        groupings[""] = sorted_parameters
    end

    # Open the TeX file
    basename = use_mode ? "prior_posterior_mode" : "prior_posterior_means"
    if !isempty(subset_string)
        basename *= "_sub=$(subset_string)"
    end
    table_out = tablespath(m, "estimate", "$basename.tex")
    if !isempty(outdir)
        outfile = replace(outfile, dirname(outfile) => outdir)
    end
    fid = open(table_out, "w")

    # Write header
    write_table_preamble(fid)
    @printf fid "\\vspace*{.5cm}\n"
    @printf fid "\\begin{longtable}{ccc}\n"
    if caption
        if use_mode
            @printf fid "\\caption{Parameter Estimates: Prior Mean and Posterior Mode}\n"
        else
            @printf fid "\\caption{Parameter Estimates: Prior and Posterior Means}\n"
        end
    end
    @printf fid "\\\\ \\hline\n"
    @printf fid "Parameter & Prior & Posterior\n"
    @printf fid "\\\\ \\hline\n"
    @printf fid "\\endhead\n"
    @printf fid "\\hline\n"
    @printf fid "\\endfoot\n"

    # Write results
    # sorted_parameters = sort(m.parameters, by = (x -> x.key))

    for group_desc in keys(groupings)
        params = groupings[group_desc]

        # Write grouping description if not empty
        if !isempty(group_desc)
            @printf fid "\\multicolumn{7}{c}{\\textit{%s}} \\\\[3pt]\n" group_desc
        end

        for param in params
            index = m.keys[param.key]

            post_value = if param.fixed
                param.value
            else
                prior = get(param.prior)
                isa(prior, RootInverseGamma) ? prior.τ : mean(prior)
            end

            @printf fid "\$ %4.99s\$ & " param.tex_label
            @printf fid "%8.3f & " post_value
            @printf fid "\\%8.3f \\\\\n" post_values[index]
        end
    end

    # Write footer
    write_table_postamble(fid)

    # Close file
    close(fid)
end

"""
```
find_density_bands(draws::Matrix, percent::AbstractFloat; minimize::Bool=true)
```

Returns a `2` x `cols(draws)` matrix `bands` such that `percent` of the mass of `draws[:,i]`
is above `bands[1,i]` and below `bands[2,i]`.

### Arguments

- `draws`: `ndraws` by `nperiods` matrix of parameter draws (from Metropolis-Hastings, for example)
- `percent`: percent of data within bands (e.g. .9 to get 90% of mass within bands)

### Optional Arguments

- `minimize`: if `true`, choose shortest interval, otherwise just chop off lowest and
  highest (percent/2)
"""
function find_density_bands(draws::AbstractArray, percent::T;
                            minimize::Bool = true) where {T<:AbstractFloat}
    if !(0 <= percent <= 1)
        error("percent must be between 0 and 1")
    end

    ndraws, nperiods = size(draws)

    if ndraws == 1
        band = repeat(draws, outer=(2, 1))
        return band
    end

    band = zeros(2, nperiods)
    n_in_band  = round(Int, percent * ndraws)  # number of draws in the band

    for i in 1:nperiods

        # Sort response for parameter i such that 1st element is largest
        draw_variable_i = draws[:,i]
        sort!(draw_variable_i, rev=true)

        # Search to find the shortest interval containing `percent` of
        # the mass `low` is the index of the largest draw in the band
        # (but the first index to take in `draw_variable_i`, `high` is
        # the smallest (but the highest index to take)

        low = if minimize

            low         = 1
            done        = 0
            j           = 2
            minwidth = draw_variable_i[1] - draw_variable_i[n_in_band]

            while j <= (ndraws - n_in_band + 1)

                newwidth = draw_variable_i[j] - draw_variable_i[j + n_in_band - 1]

                if newwidth < minwidth
                    low = j
                    minwidth = newwidth
                end

                j += 1
            end

            low
        else
            # Chop off lowest and highest percent/2
            ndraws - n_in_band - round(Int, floor(.5*(ndraws-n_in_band)))
        end

        high = low + n_in_band - 1

        band[2,i] = draw_variable_i[low]
        band[1,i] = draw_variable_i[high]
    end
    return band
end

"""
```
find_density_bands(draws::Matrix, percents::Vector{T}; minimize::Bool=true) where T<:AbstractFloat
```

Returns a `2` x `cols(draws)` matrix `bands` such that `percent` of the mass of `draws[:,i]`
is above `bands[1,i]` and below `bands[2,i]`.

### Arguments

- `draws`: Matrix of parameter draws (from Metropolis-Hastings, for example)
- `percent`: percent of data within bands (e.g. .9 to get 90% of mass within bands)

### Optional Arguments

- `minimize`: if `true`, choose shortest interval, otherwise just chop off lowest and
  highest (percent/2)
"""
function find_density_bands(draws::AbstractArray, percents::Vector{T};
                            minimize::Bool = true) where {T<:AbstractFloat}
    bands = DataFrame()

    for p in percents
        out = find_density_bands(draws, p, minimize = minimize)
        bands[!, Symbol("$(100*p)% UB")] = vec(out[2,:])
        bands[!, Symbol("$(100*p)% LB")] = vec(out[1,:])
    end
    bands
end

function write_table_preamble(fid::IOStream)
    @printf fid "\\documentclass[12pt]{article}\n"
    @printf fid "\\usepackage{booktabs}\n"
    @printf fid "\\usepackage[justification=centering]{caption}\n"
    @printf fid "\\usepackage[margin=1in]{geometry}\n"
    @printf fid "\\usepackage{longtable}\n"
    @printf fid "\\usepackage{graphicx}\n"
    @printf fid "\\usepackage{cellspace}\n"
    @printf fid "\\setlength\\cellspacetoplimit{7pt}\n"
    @printf fid "\\setlength\\cellspacebottomlimit{7pt}\n"
    @printf fid "\\begin{document}\n"
    @printf fid "\\pagestyle{empty}\n"
end

# `small`: Whether to print an additional curly bracket after "\end{longtable}" (necessary if
# the table is enclosed by "\small{}")
function write_table_postamble(fid::IOStream; small::Bool=false, tabular::Bool=false)
    if small
        if tabular
            @printf fid "\\end{tabular}}\n"
        else
            @printf fid "\\end{longtable}}\n"
        end
    else
        if tabular
            @printf fid "\\end{tabular}\n"
        else
            @printf fid "\\end{longtable}\n"
        end
    end

    @printf fid "\\end{document}"
end

"""
```
sample_λ(m, pred_dens, θs, T = -1; parallel = true) where S<:AbstractFloat
sample_λ(m, pred_dens, T = -1; parallel = true) where S<:AbstractFloat
```

Computes and samples from the conditional density p(λ_t|θ, I_t, P) for
particle in `θs`, which represents the posterior distribution. The sampled
λ particles represent the posterior distribution p(λ_{t|t} | I_t, P).

If no posterior distribution is passed in, then the function computes
the distribution of λ_{t|t} for a static pool.

### Inputs

- `m::PoolModel{S}`: `PoolModel` object
- `pred_dens::Matrix{S}`: matrix of predictive densities
- `θs::Matrix{S}`: matrix of particles representing posterior distribution of θ
- `T::Int64`: final period for tempered particle filter

where `S<:AbstractFloat`.

### Keyword Argument

- `parallel::Bool`: use parallel computing to compute and sample draws of λ

### Outputs

- `λ_sample::Vector{Float64}`: sample of draws of λs; together with (θ,λ) represents a joint density

```
"""
function sample_λ(m::PoolModel{S}, pred_dens::Matrix{S}, θs::Matrix{S}, T::Int64 = -1;
                  parallel::Bool = false,
                  tuning0::Dict{Symbol,Any} = Dict{Symbol,Any}()) where S<:AbstractFloat
    # Check size and orientation of θs is correct: assume particle_num x parameter_num
    if length(m.parameters) != size(θs,2)
        error("number of parameters in PoolModel do not match number of parameters in matrix of posterior draws of θ")
    end

    # Initialize necessary objects
    if T <= 0
         error("T must be positive") # No period provided or is invalid
    end
    tuning = isempty(tuning0) ? deepcopy(get_setting(m, :tuning)) : deepcopy(tuning0)
    tuning[:get_t_particle_dist] = true
    tuning[:parallel] = false
    tuning[:allout] = false

    # Sample from p(λ|θ, I_t^P, P) for each θ in posterior
    data = (T == 1) ? reshape(pred_dens[:,1], 2, 1) : pred_dens[:,1:T]
    if parallel
        # Send variables to workers to avoid issues with serialization
        # across workers with different Julia system images
        sendto(workers(), m = m)
        sendto(workers(), data = data)
        sendto(workers(), θs = θs)
        sendto(workers(), tuning = tuning)

        # Sample from λ distribution given θ
        λ_sample = @sync @distributed (vcat) for i in 1:size(θs,1)
            sample_λ(m, data, vec(θs[i,:]), tuning)
        end
        if sum(λ_sample) == 0
            error("Sums to zero")
        end
        λ_sample = Array(λ_sample)
    else
        # Same as above but everything is serialized
        λ_sample = Vector{Float64}(undef, size(θs,1))
        for i in 1:size(θs,1)
            λ_sample[i] = sample_λ(m, data, vec(θs[i,:]), tuning)
        end
    end

    return λ_sample
end

# This function actually does the sampling for a given θ,
# but we provide a wrapper for an easier user experience
function sample_λ(m::PoolModel{S}, data::Matrix{S}, θ::Vector{S},
                  tuning::Dict{Symbol,Any}) where S<:AbstractFloat
    update!(m, θ)
    loglik, λ_particles, λ_weights = DSGE.filter(m, data; tuning = tuning)
    return λ_particles[size(data,2)][1,DSGE.sample(DSGE.Weights(λ_weights[:,end]))]
end

function sample_λ(m::PoolModel{S}, pred_dens::Matrix{S}, T::Int64 = -1;
                  parallel::Bool = false,
                  filestring_addl::Vector{String} = Vector{String}(undef,0),
                  tuning0::Dict{Symbol,Any} = Dict{Symbol,Any}()) where S<:AbstractFloat

    # Check m is static
    if get_setting(m, :weight_type) != :static
        error("PoolModel is not static. Set the keyword argument weight_type = :static to create a static PoolModel object.")
    end

    # Initialize necessary objects
    if T <= 0
         error("T must be positive") # No period provided or is invalid
    end
    tuning = isempty(tuning0) ? deepcopy(get_setting(m, :tuning)) : deepcopy(tuning0)
    tuning[:get_t_particle_dist] = true
    tuning[:parallel] = parallel
    tuning[:allout] = false
    orig_samp_method = get_setting(m, :sampling_method)
    m <= Setting(:sampling_method, :MH)

    # Compute posterior from a static pool
    data = (T == 1) ? reshape(pred_dens[:,1], 2, 1) : pred_dens[:,1:T] # make sure it is matrix
    estimate(m, data; filestring_addl = filestring_addl, proposal_covariance = ones(1,1))
    m <= Setting(:sampling_method, orig_samp_method)

    return h5read(rawpath(m, "estimate", "mhsave.h5", filestring_addl), "mhparams")
end

"""
```
compute_Eλ(m, h, λvec, θmat = [], weights = [];
    current_period = true, parallel = true) where T<:AbstractFloat
```

Computes and samples from the conditional density p(λ_t|θ, I_t, P) for
particle in `θs`, which represents the posterior distribution.

### Inputs

- `m::PoolModel{T}`: `PoolModel` object
- `h::Int64`: forecast horizon
- `λvec::Vector{T}`: vector of particles of λ samples from (θ,λ) joint distribution
- `θmat::Matrix{T}': matrix of posterior parameter samples
- `weights::Vector{T}`: weights of λ particles, defaults to equal weights

### Keyword Argument

- `current_period::Bool`: compute Eλ for current period t
- `parallel::Bool`: use parallel computing to compute and sample λ
- `get_dpp_pred_dens::Bool`: compute predictive densities according to dynamic prediction pools

### Outputs

- `λhat_tplush::Float64`: E[λ_{t+h|t} | I_t^P, P]
- `λhat_t::Float64`: E[λ_{t|t} | I_t^P, P]
```
"""
function compute_Eλ(m::PoolModel{T}, h::Int64, λvec::Vector{T},
                    θmat::Matrix{T} = Matrix{Float64}(undef,0,0),
                    weights::Vector{T} = Vector{Float64}(undef,0);
                    current_period::Bool = true, parallel::Bool = false) where T<:AbstractFloat

    # Set up
    if isempty(weights)
        weights = ones(length(λvec)) # assume equal weights
    end
    update_param = !isempty(θmat)
    if update_param # save so we can reset the parameters of the PoolModel
        old_param = map(x -> x.value, m.parameters)
    end
    λhat_t = if current_period mean(λvec .* weights) end # compute expected lambda in current period t

    # Push forward states and compute mean
    if parallel
        # Send data to all workers
        sendto(workers(), m = m)
        sendto(workers(), h = h)
        if update_param
            sendto(workers(), θmat = θmat)

            # Propagate forward!
            λ_vec = @sync @distributed (vcat) for i in 1:length(λvec)
                propagate_λ(λvec[i], h, m, vec(θmat[i,:]))
            end
        else
            λ_vec = @sync @distributed (vcat) for i in 1:length(λvec)
                propagate_λ(λvec[i], h, m)
            end
        end
    else
        if update_param
            λ_vec = map(i -> propagate_λ(λvec[i], h, m, vec(θmat[i,:])), collect(1:length(λvec)))
        else
            λ_vec = map(x -> propagate_λ(x, h, m), λvec)
        end
    end
    λhat_tplush = mean(λ_vec .* weights)

    if update_param
        # Update m to take original parameters
        DSGE.update!(m, old_param)
    end

    if current_period
        return λhat_tplush, λhat_t
    else
        return λhat_tplush
    end
end

"""
```
propgate_λ(λvec, h, m, θvec) where T<:AbstractFloat
```

Propagates a λ particle h periods forward.

### Inputs

- `λ::T`: λ sample from (θ,λ) joint distribution
- `h::Int64`: forecast horizon
- `m::PoolModel`: PoolModel object
- `θvec::Vector{T}`: optional vector of parameters to update PoolModel

```
"""
function propagate_λ(λ::T, h::Int64, m::PoolModel,
                     θvec = Vector{Float64}(undef,0)) where T<:AbstractFloat
    if !isempty(θvec)
        update!(m, θvec)
    end
    Φ, ~, ~ = solve(m)
    for j in 1:h
        λ = Φ([λ; 1 - λ], [0.])[1]
    end
    return λ
end


text:
### A Pluto.jl notebook ###
# v0.14.0

using Markdown
using InteractiveUtils

# ╔═╡ 15e16bb0-3560-11ec-238e-e3018fb17cf9
begin

# Import packages used

using Distributed
using DelimitedFiles
using DataFrames
using Optim
using CSV
using LinearAlgebra
using Statistics
using Distributions
using MAT
using ForwardDiff
using SparseArrays

end

# ╔═╡ 9c0d53a0-38c3-11ec-0d4d-1d549264f0f0
# Replication of DellaVigna, List, Malmendier and Rao, 2017, "Voting to Tell Others"
# The code in this notebook performs the benchmark estimates with heterogeneous auxiliary parameters (Table 3, column 1)

# ╔═╡ 92cfa5a2-5ce9-11ec-0e3e-61fe154ce557
# Authors:
# Massimiliano Pozzi (Bocconi University, pozzi.massimiliano@studbocconi.it)
# Salvatore Nunnari (Bocconi University, https://snunnari.github.io, salvatore.nunnari@unibocconi.it)

# ╔═╡ b34d84a0-5ce9-11ec-1316-39f597588bef
# This notebook works with both Neptune and Pluto; the "begin" and "end" tags in each cell are necessary for Pluto notebooks but do not play any role in Neptune

# ╔═╡ 73767a00-5db9-11ec-1bef-b9e6f6e711cc
# This notebook was tested with the following packages versions:

# Pozzi: julia 1.6.1, Pluto 0.14.7, Neptune 0.14.0, Distributions 0.25.17, DataFrames 1.2.2, Optim 1.4.1, CSV 0.9.5, MAT 0.10.1, ForwardDiff 0.10.20

# ╔═╡ c6a68930-3588-11ec-08e1-7d0922b9816b
begin

# 1. Data Cleaning and Data Preparation

# We import the relevant 100 empirical moments and variance-covariance matrix. In addition to those moments we manually add the baseline turnout and its standard error. For the benchmark specification we then use 101 moments.

file = matopen(raw"../input/Moments.mat")
emp_moments = vec(read(file, "Moments"))
emp_moments = [emp_moments; 0.6000] # empirical moments. We manually add baseline                                           turnout

# We now import the variance-covariance matrix of the empirical moments. sparse to use blockdiag in the next lines

VCcontrol = sparse(read(file, "VCcontrol"))

# We need to add an elements on the diagonal of VCcontrol. Matrix to go back from sparse

VCcontrol = Matrix(blockdiag(VCcontrol,sparse(diagm([0.0109^2]))))

# W is a diagonal weighting matrix. This is the inverse of the diagonal of VCcontrol. We use 101 moments in the benchmark estimates, so this is a 101x101 diagonal matrix

W = inv(Diagonal(Diagonal(VCcontrol)))

end;

# ╔═╡ c6731d20-3588-11ec-3127-51520038c40c
# 2. Briefly define the model and the estimation strategy

# The estimation method is simulated minimum distance, where we impose simulated moments being equal to the empirical moments observed in the data. These simulated moments come from the fact that some of the parameters in the model are assumed to be stochastic and heterogeneous. What we will be minimizing is the weighted sum of squared distances between simulated and empirical moments.

# For more information about the model please refer to the python notebook or section 2 and 5 in the paper

# This function computes the relevant 101 moments by simulating N individuals.
# rand_set is a vector of vectors of random draws for the simulated parameters

function voteSimEndogenousVoting_vary(parameters, rand_set)

# Parameters:
# h0  = baseline probability of being at home
# r   = prob of seeing flyer
# eta = elasticity of response to sorting in and out of the house
# s   = how much you like doing a survey
# S   = social pressure cost of doing a survey
# sv  = value of saying you voted (appearing as a voter)
# sn  = value of saying you didn't vote (appearing as a non-voter)
# rho = correlation between sv and sn (assumed to be 0 in the benchmark estimates)
# eps = all other reasons to vote

N=5.4    # Number of times asked if you voted in the Congressional elections
N_P=10.1 # Times asked in presidential

# Auxiliary parameters. These parameters vary between voters and non-voters. sub v for voters and sub nv for non-voters

h0_v	=	parameters[1]
h0_nv	=	parameters[2]
r_v	    =	parameters[3]
r_nv	=	parameters[4]
eta_v	=	parameters[5]
eta_nv	=	parameters[6]
mu_s_v	=	parameters[7]
mu_s_nv	=	parameters[8]
sigma_s_v	=	parameters[9]
sigma_s_nv	=	parameters[10]
S_svy_v	    =	parameters[11]
S_svy_nv	=	parameters[12]
timeval_v	=	parameters[13]
timeval_nv	=	parameters[14]

# Other parameters. Here they don't vary between voters and non-voters
mu_sv	=	parameters[15]
mu_sn	=	parameters[16]
sigma_svn = parameters[17]  # sigma_sv = sigma_sn in benchmark estimates
L = parameters[18]
mu_eps = parameters[19]
sigma_eps = parameters[20]
rho	= 0                     # Set to zero in benchmark

# Draws from a standard normal. The dimension of these vectors is equal to the number of simulated individuals

rand_set1 = rand_set[1] # used to compute s
rand_set2 = rand_set[2] # used to compute s_v
rand_set3 = rand_set[3] # used to compute s_n
rand_set4 = rand_set[4] # used to compute eps

# Simulate epsilons ("other" reasons to vote)

eps = mu_eps .+ sigma_eps.*rand_set2

# Simulate sv and sn

sv = mu_sv .+ sigma_svn.*rand_set3
sn = mu_sn .+ sigma_svn.*rand_set4

# !!!Vote or not!!! Does our simulated individual vote or not?

#Net expected utility of voting (relative to not voting)

sigVal = (max.(sv,sn.-L) .- max.(sn,sv.-L))# net sig utility voter - non-voter (1 ask)
sigVal_x_N = N.*sigVal                     # asked N times
utilityVoting = sigVal_x_N .+ eps          # net utility
voted = utilityVoting .>0                  # Vote if net utility > 0

# Return a list of indices (e.g., 1, 3, 4) identifying voters and non-voters

voterIndex    = findall(==(1), voted)
nonvoterIndex = findall(==(0), voted)

# Share who turn out in the control group (no GOTV intervention)

Turnout_control = mean(voted)

# Make vectors with voter and non-voter parameters based on whether voted in the control experiment

h0 = voted.*h0_v .+ (1 .-voted).*h0_nv
r = voted.*r_v .+ (1 .-voted).*r_nv
eta = voted.*eta_v .+ (1 .-voted).*eta_nv
S_svy = voted.*S_svy_v .+ (1 .-voted).*S_svy_nv
timeval = voted.*timeval_v .+ (1 .-voted).*timeval_nv

# GOTV (use parameters as assigned by turnout in the control). Assumes intervention = expect to be asked h0 times more

# net utility of voting if seen the GOTV flyer
sigVal_x_N_GOTV = (N.+h0).*sigVal
utilityVoting_GOTV = sigVal_x_N_GOTV .+ eps
voted_GOTV = utilityVoting_GOTV .>0

# Share who turn out with the GOTV intervention assume everyone sees the flyer, counts as N+h0 times asked
Turnout_GOTV = mean(voted_GOTV)

# PRESIDENTIAL
sigVal_x_N_P = N_P.*sigVal; # asked N_P times
sigVal_x_N_P_GOTV = (N_P.+h0).*sigVal # asked N_P times

utilityVoting_P = sigVal_x_N_P .+ eps
utilityVoting_P_GOTV = sigVal_x_N_P_GOTV .+ eps

voted_P = utilityVoting_P .>0
voted_P_GOTV = utilityVoting_P_GOTV .>0

# Turnout
Turnout_P_control = mean(voted_P)
Turnout_P_GOTV = mean(voted_P_GOTV)

# Draw random s, Utilities of doing the survey for voters/non-voters

# Simulate utility of doing a 0d10m survey
s = voted.*mu_s_v .+ voted.*sigma_s_v.*rand_set1 .+ (1 .-voted).*mu_s_nv .+ (1 .-voted).*sigma_s_nv.*rand_set1

# Values of survey incentives (relative to 0d10m) XdYm = X dollars and Y min
D_0d10m_v = 0
D_0d5m_v = timeval_v*5/60
D_10d5m_v = 10+timeval_v*5/60
D_10d10m_v = 10

D_0d10m_nv = 0
D_0d5m_nv = timeval_nv*5/60
D_10d5m_nv = 10+timeval_nv*5/60
D_10d10m_nv = 10

# Extra incentive if say "not vote". 5m survey: +1m + $5. 10m survey: -8m
I_5d1m_v = 5 .-timeval_v*1/60
I_8m_v = timeval_v*8/60

I_5d1m_nv = 5 .-timeval_nv*1/60
I_8m_nv = timeval_nv*8/60

# Lying if asked

# UtilVotingQuestion = utility you get from being asked one time (max of lie or not lie)
wouldLieIfAsked = voted.*(sn.-L .>sv) .+ (1 .-voted).*(sv.-L .>sn)
utilVotingQuestion = voted.*max.(sn.-L,sv) .+ (1 .-voted).*max.(sv.-L,sn)

# Response to incentives to say "not vote"
wouldLieIfAsked_5d1m = voted.*(sn.-L.+I_5d1m_v .>sv) .+ (1 .-voted).*(sv.-L.>sn.+I_5d1m_nv)
wouldLieIfAsked_8m = voted.*(sn.-L.+I_8m_v .>sv) .+ (1 .-voted).*(sv.-L .>sn.+I_8m_nv)

# !!!!! Compute Moments !!!!!

# Utility from doing survey

# NF= no flyer, F = flyer, FV= flyer + voting, OO = opt-out, OOV = opt-out + voting

# Anticipated utility from doing survey and voting survey (VF or I)
util_svyOnly_0d5m = s .+ voted.*D_0d5m_v .+ (1 .-voted).*D_0d5m_nv
util_svyPlusVotingQues_0d5m = util_svyOnly_0d5m .+ utilVotingQuestion

util_svyOnly_10d10m = s .+ voted.*D_10d10m_v .+ (1 .-voted).*D_10d10m_nv
util_svyPlusVotingQues_10d10m = util_svyOnly_10d10m .+ utilVotingQuestion

util_svyOnly_10d5m = s .+ voted.*D_10d5m_v .+ (1 .-voted).*D_10d5m_nv
util_svyPlusVotingQues_10d5m = util_svyOnly_10d5m .+ utilVotingQuestion

# If asked, do survey if greater than the social pressure cost. NI = not informed that survey is about voting, I=informed (VF or I)
doesSvyIfAsked_NI_0d5m = util_svyOnly_0d5m .> -S_svy
doesSvyIfAsked_I_0d5m = util_svyPlusVotingQues_0d5m .> -S_svy

doesSvyIfAsked_NI_10d10m = util_svyOnly_10d10m .> -S_svy
doesSvyIfAsked_I_10d10m = util_svyPlusVotingQues_10d10m .> -S_svy

doesSvyIfAsked_NI_10d5m = util_svyOnly_10d5m .> -S_svy
doesSvyIfAsked_I_10d5m = util_svyPlusVotingQues_10d5m .> -S_svy

# Anticipated utility given you are asked to do the survey
anticipatedUtil_Svy_NI_0d5m = max.(util_svyOnly_0d5m,-S_svy)
anticipatedUtil_Svy_I_0d5m = max.(util_svyPlusVotingQues_0d5m,-S_svy)

anticipatedUtil_Svy_NI_10d10m = max.(util_svyOnly_10d10m,-S_svy)
anticipatedUtil_Svy_I_10d10m = max.(util_svyPlusVotingQues_10d10m,-S_svy)

anticipatedUtil_Svy_NI_10d5m = max.(util_svyOnly_10d5m,-S_svy)
anticipatedUtil_Svy_I_10d5m = max.(util_svyPlusVotingQues_10d5m,-S_svy)

# Opt-out if anticipated utility is negative
optsOutIfSees_OO_0d5m  = anticipatedUtil_Svy_NI_0d5m .< 0
optsOutIfSees_OOV_0d5m = anticipatedUtil_Svy_I_0d5m .< 0

optsOutIfSees_OO_10d10m  = anticipatedUtil_Svy_NI_10d10m .< 0
optsOutIfSees_OOV_10d10m = anticipatedUtil_Svy_I_10d10m .< 0

optsOutIfSees_OO_10d5m  = anticipatedUtil_Svy_NI_10d5m .< 0
optsOutIfSees_OOV_10d5m = anticipatedUtil_Svy_I_10d5m .< 0

# Choosing probability of being at home is bounded between 0 and 1

hStar_F_0d5m = max.(0,min.(1,h0.+eta.*anticipatedUtil_Svy_NI_0d5m))
hStar_FV_0d5m = max.(0,min.(1,h0.+eta.*anticipatedUtil_Svy_I_0d5m))

hStar_F_10d10m = max.(0,min.(1,h0.+eta.*anticipatedUtil_Svy_NI_10d10m))
hStar_FV_10d10m = max.(0,min.(1,h0.+eta.*anticipatedUtil_Svy_I_10d10m))

hStar_F_10d5m = max.(0,min.(1,h0+eta.*anticipatedUtil_Svy_NI_10d5m))
hStar_FV_10d5m = max.(0,min.(1,h0+eta.*anticipatedUtil_Svy_I_10d5m))

# Separate Voters and Nonvoters

# Split into separate vectors of voters and non-voter vectors note: they will be of different length

# Voters
hStar_F_0d5m_v = hStar_F_0d5m[voterIndex]
hStar_FV_0d5m_v = hStar_FV_0d5m[voterIndex]
doesSvyIfAsked_NI_0d5m_v = doesSvyIfAsked_NI_0d5m[voterIndex]
doesSvyIfAsked_I_0d5m_v = doesSvyIfAsked_I_0d5m[voterIndex]
optsOutIfSees_OO_0d5m_v=optsOutIfSees_OO_0d5m[voterIndex]
optsOutIfSees_OOV_0d5m_v=optsOutIfSees_OOV_0d5m[voterIndex]

hStar_F_10d10m_v = hStar_F_10d10m[voterIndex]
hStar_FV_10d10m_v = hStar_FV_10d10m[voterIndex]
doesSvyIfAsked_NI_10d10m_v = doesSvyIfAsked_NI_10d10m[voterIndex]
doesSvyIfAsked_I_10d10m_v = doesSvyIfAsked_I_10d10m[voterIndex]
optsOutIfSees_OO_10d10m_v=optsOutIfSees_OO_10d10m[voterIndex]
optsOutIfSees_OOV_10d10m_v=optsOutIfSees_OOV_10d10m[voterIndex]

hStar_F_10d5m_v = hStar_F_10d5m[voterIndex]
hStar_FV_10d5m_v = hStar_FV_10d5m[voterIndex]
doesSvyIfAsked_NI_10d5m_v = doesSvyIfAsked_NI_10d5m[voterIndex]
doesSvyIfAsked_I_10d5m_v = doesSvyIfAsked_I_10d5m[voterIndex]
optsOutIfSees_OO_10d5m_v=optsOutIfSees_OO_10d5m[voterIndex]
optsOutIfSees_OOV_10d5m_v=optsOutIfSees_OOV_10d5m[voterIndex]

wouldLieIfAsked_v = wouldLieIfAsked[voterIndex]
wouldLieIfAsked_5d1m_v = wouldLieIfAsked_5d1m[voterIndex]
wouldLieIfAsked_8m_v=wouldLieIfAsked_8m[voterIndex]

# Non-voters
hStar_F_0d5m_nv = hStar_F_0d5m[nonvoterIndex]
hStar_FV_0d5m_nv = hStar_FV_0d5m[nonvoterIndex]
doesSvyIfAsked_NI_0d5m_nv = doesSvyIfAsked_NI_0d5m[nonvoterIndex]
doesSvyIfAsked_I_0d5m_nv = doesSvyIfAsked_I_0d5m[nonvoterIndex]
optsOutIfSees_OO_0d5m_nv=optsOutIfSees_OO_0d5m[nonvoterIndex]
optsOutIfSees_OOV_0d5m_nv=optsOutIfSees_OOV_0d5m[nonvoterIndex]

hStar_F_10d10m_nv = hStar_F_10d10m[nonvoterIndex]
hStar_FV_10d10m_nv = hStar_FV_10d10m[nonvoterIndex]
doesSvyIfAsked_NI_10d10m_nv = doesSvyIfAsked_NI_10d10m[nonvoterIndex]
doesSvyIfAsked_I_10d10m_nv = doesSvyIfAsked_I_10d10m[nonvoterIndex]
optsOutIfSees_OO_10d10m_nv=optsOutIfSees_OO_10d10m[nonvoterIndex]
optsOutIfSees_OOV_10d10m_nv=optsOutIfSees_OOV_10d10m[nonvoterIndex]

hStar_F_10d5m_nv = hStar_F_10d5m[nonvoterIndex]
hStar_FV_10d5m_nv = hStar_FV_10d5m[nonvoterIndex]
doesSvyIfAsked_NI_10d5m_nv = doesSvyIfAsked_NI_10d5m[nonvoterIndex]
doesSvyIfAsked_I_10d5m_nv = doesSvyIfAsked_I_10d5m[nonvoterIndex]
optsOutIfSees_OO_10d5m_nv=optsOutIfSees_OO_10d5m[nonvoterIndex]
optsOutIfSees_OOV_10d5m_nv=optsOutIfSees_OOV_10d5m[nonvoterIndex]

wouldLieIfAsked_nv = wouldLieIfAsked[nonvoterIndex]
wouldLieIfAsked_5d1m_nv = wouldLieIfAsked_5d1m[nonvoterIndex]
wouldLieIfAsked_8m_nv=wouldLieIfAsked_8m[nonvoterIndex]

# !!!!! Disaggregated Moments !!!!!

# Voters

# PH = probability of being at home
PH_NF_0d5m_v = h0_v;
PH_F_0d5m_v = (1 .-r_v)*h0_v .+ r_v.*mean(hStar_F_0d5m_v);
PH_FV_0d5m_v = (1 .-r_v)*h0_v .+ r_v.*mean(hStar_FV_0d5m_v);
PH_OO_0d5m_v = (1 .-r_v)*h0_v .+ r_v.*mean((1 .-optsOutIfSees_OO_0d5m_v).*hStar_F_0d5m_v);
PH_OOV_0d5m_v = (1 .-r_v)*h0_v .+ r_v.*mean((1 .-optsOutIfSees_OOV_0d5m_v).*hStar_FV_0d5m_v);

PH_NF_10d10m_v = h0_v;
PH_F_10d10m_v = (1 .-r_v)*h0_v .+ r_v.*mean(hStar_F_10d10m_v);
PH_FV_10d10m_v = (1 .-r_v)*h0_v .+ r_v.*mean(hStar_FV_10d10m_v);
PH_OO_10d10m_v = (1 .-r_v)*h0_v .+ r_v.*mean((1 .-optsOutIfSees_OO_10d10m_v).*hStar_F_10d10m_v);
PH_OOV_10d10m_v = (1 .-r_v)*h0_v .+ r_v.*mean((1 .-optsOutIfSees_OOV_10d10m_v).*hStar_FV_10d10m_v);

PH_NF_10d5m_v = h0_v;
PH_F_10d5m_v = (1 .-r_v)*h0_v .+ r_v.*mean(hStar_F_10d5m_v);
PH_FV_10d5m_v = (1 .-r_v)*h0_v .+ r_v.*mean(hStar_FV_10d5m_v);
PH_OO_10d5m_v = (1 .-r_v)*h0_v .+ r_v.*mean((1 .-optsOutIfSees_OO_10d5m_v).*hStar_F_10d5m_v);
PH_OOV_10d5m_v = (1 .-r_v)*h0_v .+ r_v.*mean((1 .-optsOutIfSees_OOV_10d5m_v).*hStar_FV_10d5m_v);

# PSV=unconditional prob of doing the survey (not cond on opening door). PSV < PH mechanically

# 0d5m
PSV_NF_NI_0d5m_v = h0_v*mean(doesSvyIfAsked_NI_0d5m_v);
PSV_NF_I_0d5m_v = h0_v*mean(doesSvyIfAsked_I_0d5m_v);

PSV_F_NI_0d5m_v = (1 .-r_v).*PSV_NF_NI_0d5m_v .+ r_v.*mean(hStar_F_0d5m_v.*doesSvyIfAsked_NI_0d5m_v);
PSV_F_I_0d5m_v = (1 .-r_v)*PSV_NF_I_0d5m_v .+ r_v.*mean(hStar_F_0d5m_v.*doesSvyIfAsked_I_0d5m_v);

PSV_FV_NI_0d5m_v = (1 .-r_v)*PSV_NF_NI_0d5m_v .+ r_v.*mean(hStar_FV_0d5m_v.*doesSvyIfAsked_I_0d5m_v);
PSV_FV_I_0d5m_v = (1 .-r_v)*PSV_NF_I_0d5m_v .+ r_v.*mean(hStar_FV_0d5m_v.*doesSvyIfAsked_I_0d5m_v);

PSV_OO_NI_0d5m_v = (1 .-r_v)*PSV_NF_NI_0d5m_v .+ r_v*mean((1 .-optsOutIfSees_OO_0d5m_v).*hStar_F_0d5m_v.*doesSvyIfAsked_NI_0d5m_v);
PSV_OO_I_0d5m_v = (1 .-r_v)*PSV_NF_I_0d5m_v .+ r_v*mean((1 .-optsOutIfSees_OO_0d5m_v).*hStar_F_0d5m_v.*doesSvyIfAsked_I_0d5m_v);

PSV_OOV_NI_0d5m_v = (1 .-r_v)*PSV_NF_NI_0d5m_v .+ r_v*mean((1 .-optsOutIfSees_OOV_0d5m_v).*hStar_FV_0d5m_v.*doesSvyIfAsked_I_0d5m_v);
PSV_OOV_I_0d5m_v = (1 .-r_v)*PSV_NF_I_0d5m_v .+ r_v*mean((1 .-optsOutIfSees_OOV_0d5m_v).*hStar_FV_0d5m_v.*doesSvyIfAsked_I_0d5m_v);

# 10d10m
PSV_NF_NI_10d10m_v = h0_v*mean(doesSvyIfAsked_NI_10d10m_v);
PSV_NF_I_10d10m_v = h0_v*mean(doesSvyIfAsked_I_10d10m_v);

PSV_F_NI_10d10m_v = (1-r_v)*PSV_NF_NI_10d10m_v + r_v*mean(hStar_F_10d10m_v.*doesSvyIfAsked_NI_10d10m_v);
PSV_F_I_10d10m_v = (1-r_v)*PSV_NF_I_10d10m_v + r_v*mean(hStar_F_10d10m_v.*doesSvyIfAsked_I_10d10m_v);

PSV_FV_NI_10d10m_v = (1-r_v)*PSV_NF_NI_10d10m_v + r_v*mean(hStar_FV_10d10m_v.*doesSvyIfAsked_I_10d10m_v);
PSV_FV_I_10d10m_v = (1-r_v)*PSV_NF_I_10d10m_v + r_v*mean(hStar_FV_10d10m_v.*doesSvyIfAsked_I_10d10m_v);

PSV_OO_NI_10d10m_v = (1-r_v)*PSV_NF_NI_10d10m_v + r_v*mean((1 .-optsOutIfSees_OO_10d10m_v).*hStar_F_10d10m_v.*doesSvyIfAsked_NI_10d10m_v);
PSV_OO_I_10d10m_v = (1-r_v)*PSV_NF_I_10d10m_v + r_v*mean((1 .-optsOutIfSees_OO_10d10m_v).*hStar_F_10d10m_v.*doesSvyIfAsked_I_10d10m_v);

PSV_OOV_NI_10d10m_v = (1-r_v)*PSV_NF_NI_10d10m_v + r_v*mean((1 .-optsOutIfSees_OOV_10d10m_v).*hStar_FV_10d10m_v.*doesSvyIfAsked_I_10d10m_v);
PSV_OOV_I_10d10m_v = (1-r_v)*PSV_NF_I_10d10m_v + r_v*mean((1 .-optsOutIfSees_OOV_10d10m_v).*hStar_FV_10d10m_v.*doesSvyIfAsked_I_10d10m_v);

# 10d5m
PSV_NF_NI_10d5m_v = h0_v*mean(doesSvyIfAsked_NI_10d5m_v);
PSV_NF_I_10d5m_v = h0_v*mean(doesSvyIfAsked_I_10d5m_v);

PSV_F_NI_10d5m_v = (1-r_v)*PSV_NF_NI_10d5m_v + r_v*mean(hStar_F_10d5m_v.*doesSvyIfAsked_NI_10d5m_v);
PSV_F_I_10d5m_v = (1-r_v)*PSV_NF_I_10d5m_v + r_v*mean(hStar_F_10d5m_v.*doesSvyIfAsked_I_10d5m_v);

PSV_FV_NI_10d5m_v = (1-r_v)*PSV_NF_NI_10d5m_v + r_v*mean(hStar_FV_10d5m_v.*doesSvyIfAsked_I_10d5m_v);
PSV_FV_I_10d5m_v = (1-r_v)*PSV_NF_I_10d5m_v + r_v*mean(hStar_FV_10d5m_v.*doesSvyIfAsked_I_10d5m_v);

PSV_OO_NI_10d5m_v = (1-r_v)*PSV_NF_NI_10d5m_v + r_v*mean((1 .-optsOutIfSees_OO_10d5m_v).*hStar_F_10d5m_v.*doesSvyIfAsked_NI_10d5m_v);
PSV_OO_I_10d5m_v = (1-r_v)*PSV_NF_I_10d5m_v + r_v*mean((1 .-optsOutIfSees_OO_10d5m_v).*hStar_F_10d5m_v.*doesSvyIfAsked_I_10d5m_v);

PSV_OOV_NI_10d5m_v = (1-r_v)*PSV_NF_NI_10d5m_v + r_v*mean((1 .-optsOutIfSees_OOV_10d5m_v).*hStar_FV_10d5m_v.*doesSvyIfAsked_I_10d5m_v);
PSV_OOV_I_10d5m_v = (1-r_v)*PSV_NF_I_10d5m_v + r_v*mean((1 .-optsOutIfSees_OOV_10d5m_v).*hStar_FV_10d5m_v.*doesSvyIfAsked_I_10d5m_v);

# POO=prob of opting out (not conditional on seeing flyer). Scaled by baseline likelihood of being at home
POO_OO_0d5m_v =  h0_v*r_v*mean(optsOutIfSees_OO_0d5m_v);
POO_OOV_0d5m_v = h0_v*r_v*mean(optsOutIfSees_OOV_0d5m_v);

POO_OO_10d10m_v =  h0_v*r_v*mean(optsOutIfSees_OO_10d10m_v);
POO_OOV_10d10m_v = h0_v*r_v*mean(optsOutIfSees_OOV_10d10m_v);

POO_OO_10d5m_v =  h0_v*r_v*mean(optsOutIfSees_OO_10d5m_v);
POO_OOV_10d5m_v = h0_v*r_v*mean(optsOutIfSees_OOV_10d5m_v);

# Empirical moments are total lying in treatments / total doing survey in treatments

# PSVL = unconditional percent who do survey and lie. No flyer treatment only, simplifies later code

# PL=cond on agreeing to do the survey, did you lie? Incentive to lie is a surprise later (doesn't affect PH or PSV)

# 0d5m, 5d1m incentive
PSVL_NF_NI_0d5m_v = mean(h0_v.*doesSvyIfAsked_NI_0d5m_v.*wouldLieIfAsked_v);
PSVL_NF_I_0d5m_v = mean(h0_v.*doesSvyIfAsked_I_0d5m_v.*wouldLieIfAsked_v);
PSVL_NF_NI_0d5m_5d1m_v = mean(h0_v.*doesSvyIfAsked_NI_0d5m_v.*wouldLieIfAsked_5d1m_v);
PSVL_NF_I_0d5m_5d1m_v = mean(h0_v.*doesSvyIfAsked_I_0d5m_v.*wouldLieIfAsked_5d1m_v);

PSVL_F_NI_0d5m_v = (1-r_v)*PSVL_NF_NI_0d5m_v + r_v*mean(hStar_F_0d5m_v.*doesSvyIfAsked_NI_0d5m_v.*wouldLieIfAsked_v);
PSVL_F_I_0d5m_v = (1-r_v)*PSVL_NF_I_0d5m_v + r_v*mean(hStar_F_0d5m_v.*doesSvyIfAsked_I_0d5m_v.*wouldLieIfAsked_v);
PSVL_F_NI_0d5m_5d1m_v = (1-r_v)*PSVL_NF_NI_0d5m_5d1m_v + r_v*mean(hStar_F_0d5m_v.*doesSvyIfAsked_NI_0d5m_v.*wouldLieIfAsked_5d1m_v);
PSVL_F_I_0d5m_5d1m_v = (1-r_v)*PSVL_NF_I_0d5m_5d1m_v + r_v*mean(hStar_F_0d5m_v.*doesSvyIfAsked_I_0d5m_v.*wouldLieIfAsked_5d1m_v);

PSVL_FV_NI_0d5m_v = (1-r_v)*PSVL_NF_NI_0d5m_v + r_v*mean(hStar_FV_0d5m_v.*doesSvyIfAsked_I_0d5m_v.*wouldLieIfAsked_v);
PSVL_FV_I_0d5m_v = (1-r_v)*PSVL_NF_I_0d5m_v + r_v*mean(hStar_FV_0d5m_v.*doesSvyIfAsked_I_0d5m_v.*wouldLieIfAsked_v);
PSVL_FV_NI_0d5m_5d1m_v = (1-r_v)*PSVL_NF_NI_0d5m_5d1m_v + r_v*mean(hStar_FV_0d5m_v.*doesSvyIfAsked_I_0d5m_v.*wouldLieIfAsked_5d1m_v);
PSVL_FV_I_0d5m_5d1m_v = (1-r_v)*PSVL_NF_I_0d5m_5d1m_v + r_v*mean(hStar_FV_0d5m_v.*doesSvyIfAsked_I_0d5m_v.*wouldLieIfAsked_5d1m_v);

PSVL_OO_NI_0d5m_v = (1-r_v)*PSVL_NF_NI_0d5m_v + r_v*mean((1 .-optsOutIfSees_OO_0d5m_v).*hStar_F_0d5m_v.*doesSvyIfAsked_NI_0d5m_v.*wouldLieIfAsked_v);
PSVL_OO_I_0d5m_v = (1-r_v)*PSVL_NF_I_0d5m_v + r_v*mean((1 .-optsOutIfSees_OO_0d5m_v).*hStar_F_0d5m_v.*doesSvyIfAsked_I_0d5m_v.*wouldLieIfAsked_v);
PSVL_OO_NI_0d5m_5d1m_v = (1-r_v)*PSVL_NF_NI_0d5m_5d1m_v + r_v*mean((1 .-optsOutIfSees_OO_0d5m_v).*hStar_F_0d5m_v.*doesSvyIfAsked_NI_0d5m_v.*wouldLieIfAsked_5d1m_v);
PSVL_OO_I_0d5m_5d1m_v = (1-r_v)*PSVL_NF_I_0d5m_5d1m_v + r_v*mean((1 .-optsOutIfSees_OO_0d5m_v).*hStar_F_0d5m_v.*doesSvyIfAsked_I_0d5m_v.*wouldLieIfAsked_5d1m_v);

PSVL_OOV_NI_0d5m_v = (1-r_v)*PSVL_NF_NI_0d5m_v + r_v*mean((1 .-optsOutIfSees_OOV_0d5m_v).*hStar_FV_0d5m_v.*doesSvyIfAsked_I_0d5m_v.*wouldLieIfAsked_v);
PSVL_OOV_I_0d5m_v = (1-r_v)*PSVL_NF_I_0d5m_v + r_v*mean((1 .-optsOutIfSees_OOV_0d5m_v).*hStar_FV_0d5m_v.*doesSvyIfAsked_I_0d5m_v.*wouldLieIfAsked_v);
PSVL_OOV_NI_0d5m_5d1m_v = (1-r_v)*PSVL_NF_NI_0d5m_5d1m_v + r_v*mean((1 .-optsOutIfSees_OOV_0d5m_v).*hStar_FV_0d5m_v.*doesSvyIfAsked_I_0d5m_v.*wouldLieIfAsked_5d1m_v);
PSVL_OOV_I_0d5m_5d1m_v = (1-r_v)*PSVL_NF_I_0d5m_5d1m_v + r_v*mean((1 .-optsOutIfSees_OOV_0d5m_v).*hStar_FV_0d5m_v.*doesSvyIfAsked_I_0d5m_v.*wouldLieIfAsked_5d1m_v);

# 10d10m, 8m incentive
PSVL_NF_NI_10d10m_v = mean(h0_v.*doesSvyIfAsked_NI_10d10m_v.*wouldLieIfAsked_v);
PSVL_NF_I_10d10m_v = mean(h0_v.*doesSvyIfAsked_I_10d10m_v.*wouldLieIfAsked_v);
PSVL_NF_NI_10d10m_8m_v = mean(h0_v.*doesSvyIfAsked_NI_10d10m_v.*wouldLieIfAsked_8m_v);
PSVL_NF_I_10d10m_8m_v = mean(h0_v.*doesSvyIfAsked_I_10d10m_v.*wouldLieIfAsked_8m_v);

PSVL_F_NI_10d10m_v = (1-r_v)*PSVL_NF_NI_10d10m_v + r_v*mean(hStar_F_10d10m_v.*doesSvyIfAsked_NI_10d10m_v.*wouldLieIfAsked_v);
PSVL_F_I_10d10m_v = (1-r_v)*PSVL_NF_I_10d10m_v + r_v*mean(hStar_F_10d10m_v.*doesSvyIfAsked_I_10d10m_v.*wouldLieIfAsked_v);
PSVL_F_NI_10d10m_8m_v = (1-r_v)*PSVL_NF_NI_10d10m_8m_v + r_v*mean(hStar_F_10d10m_v.*doesSvyIfAsked_NI_10d10m_v.*wouldLieIfAsked_8m_v);
PSVL_F_I_10d10m_8m_v = (1-r_v)*PSVL_NF_I_10d10m_8m_v + r_v*mean(hStar_F_10d10m_v.*doesSvyIfAsked_I_10d10m_v.*wouldLieIfAsked_8m_v);

PSVL_FV_NI_10d10m_v = (1-r_v)*PSVL_NF_NI_10d10m_v + r_v*mean(hStar_FV_10d10m_v.*doesSvyIfAsked_I_10d10m_v.*wouldLieIfAsked_v);
PSVL_FV_I_10d10m_v = (1-r_v)*PSVL_NF_I_10d10m_v + r_v*mean(hStar_FV_10d10m_v.*doesSvyIfAsked_I_10d10m_v.*wouldLieIfAsked_v);
PSVL_FV_NI_10d10m_8m_v = (1-r_v)*PSVL_NF_NI_10d10m_8m_v + r_v*mean(hStar_FV_10d10m_v.*doesSvyIfAsked_I_10d10m_v.*wouldLieIfAsked_8m_v);
PSVL_FV_I_10d10m_8m_v = (1-r_v)*PSVL_NF_I_10d10m_8m_v + r_v*mean(hStar_FV_10d10m_v.*doesSvyIfAsked_I_10d10m_v.*wouldLieIfAsked_8m_v);

PSVL_OO_NI_10d10m_v = (1-r_v)*PSVL_NF_NI_10d10m_v + r_v*mean((1 .-optsOutIfSees_OO_10d10m_v).*hStar_F_10d10m_v.*doesSvyIfAsked_NI_10d10m_v.*wouldLieIfAsked_v);
PSVL_OO_I_10d10m_v = (1-r_v)*PSVL_NF_I_10d10m_v + r_v*mean((1 .-optsOutIfSees_OO_10d10m_v).*hStar_F_10d10m_v.*doesSvyIfAsked_I_10d10m_v.*wouldLieIfAsked_v);
PSVL_OO_NI_10d10m_8m_v = (1-r_v)*PSVL_NF_NI_10d10m_8m_v + r_v*mean((1 .-optsOutIfSees_OO_10d10m_v).*hStar_F_10d10m_v.*doesSvyIfAsked_NI_10d10m_v.*wouldLieIfAsked_8m_v);
PSVL_OO_I_10d10m_8m_v = (1-r_v)*PSVL_NF_I_10d10m_8m_v + r_v*mean((1 .-optsOutIfSees_OO_10d10m_v).*hStar_F_10d10m_v.*doesSvyIfAsked_I_10d10m_v.*wouldLieIfAsked_8m_v);

PSVL_OOV_NI_10d10m_v = (1-r_v)*PSVL_NF_NI_10d10m_v + r_v*mean((1 .-optsOutIfSees_OOV_10d10m_v).*hStar_FV_10d10m_v.*doesSvyIfAsked_I_10d10m_v.*wouldLieIfAsked_v);
PSVL_OOV_I_10d10m_v = (1-r_v)*PSVL_NF_I_10d10m_v + r_v*mean((1 .-optsOutIfSees_OOV_10d10m_v).*hStar_FV_10d10m_v.*doesSvyIfAsked_I_10d10m_v.*wouldLieIfAsked_v);
PSVL_OOV_NI_10d10m_8m_v = (1-r_v)*PSVL_NF_NI_10d10m_8m_v + r_v*mean((1 .-optsOutIfSees_OOV_10d10m_v).*hStar_FV_10d10m_v.*doesSvyIfAsked_I_10d10m_v.*wouldLieIfAsked_8m_v);
PSVL_OOV_I_10d10m_8m_v = (1-r_v)*PSVL_NF_I_10d10m_8m_v + r_v*mean((1 .-optsOutIfSees_OOV_10d10m_v).*hStar_FV_10d10m_v.*doesSvyIfAsked_I_10d10m_v.*wouldLieIfAsked_8m_v);

# 10d5m, 5d1m incentive
PSVL_NF_NI_10d5m_v = mean(h0_v.*doesSvyIfAsked_NI_10d5m_v.*wouldLieIfAsked_v);
PSVL_NF_I_10d5m_v = mean(h0_v.*doesSvyIfAsked_I_10d5m_v.*wouldLieIfAsked_v);
PSVL_NF_NI_10d5m_5d1m_v = mean(h0_v.*doesSvyIfAsked_NI_10d5m_v.*wouldLieIfAsked_5d1m_v);
PSVL_NF_I_10d5m_5d1m_v = mean(h0_v.*doesSvyIfAsked_I_10d5m_v.*wouldLieIfAsked_5d1m_v);

PSVL_F_NI_10d5m_v = (1-r_v)*PSVL_NF_NI_10d5m_v + r_v*mean(hStar_F_10d5m_v.*doesSvyIfAsked_NI_10d5m_v.*wouldLieIfAsked_v);
PSVL_F_I_10d5m_v = (1-r_v)*PSVL_NF_I_10d5m_v + r_v*mean(hStar_F_10d5m_v.*doesSvyIfAsked_I_10d5m_v.*wouldLieIfAsked_v);
PSVL_F_NI_10d5m_5d1m_v = (1-r_v)*PSVL_NF_NI_10d5m_5d1m_v + r_v*mean(hStar_F_10d5m_v.*doesSvyIfAsked_NI_10d5m_v.*wouldLieIfAsked_5d1m_v);
PSVL_F_I_10d5m_5d1m_v = (1-r_v)*PSVL_NF_I_10d5m_5d1m_v + r_v*mean(hStar_F_10d5m_v.*doesSvyIfAsked_I_10d5m_v.*wouldLieIfAsked_5d1m_v);

PSVL_FV_NI_10d5m_v = (1-r_v)*PSVL_NF_NI_10d5m_v + r_v*mean(hStar_FV_10d5m_v.*doesSvyIfAsked_I_10d5m_v.*wouldLieIfAsked_v);
PSVL_FV_I_10d5m_v = (1-r_v)*PSVL_NF_I_10d5m_v + r_v*mean(hStar_FV_10d5m_v.*doesSvyIfAsked_I_10d5m_v.*wouldLieIfAsked_v);
PSVL_FV_NI_10d5m_5d1m_v = (1-r_v)*PSVL_NF_NI_10d5m_5d1m_v + r_v*mean(hStar_FV_10d5m_v.*doesSvyIfAsked_I_10d5m_v.*wouldLieIfAsked_5d1m_v);
PSVL_FV_I_10d5m_5d1m_v = (1-r_v)*PSVL_NF_I_10d5m_5d1m_v + r_v*mean(hStar_FV_10d5m_v.*doesSvyIfAsked_I_10d5m_v.*wouldLieIfAsked_5d1m_v);

PSVL_OO_NI_10d5m_v = (1-r_v)*PSVL_NF_NI_10d5m_v + r_v*mean((1 .-optsOutIfSees_OO_10d5m_v).*hStar_F_10d5m_v.*doesSvyIfAsked_NI_10d5m_v.*wouldLieIfAsked_v);
PSVL_OO_I_10d5m_v = (1-r_v)*PSVL_NF_I_10d5m_v + r_v*mean((1 .-optsOutIfSees_OO_10d5m_v).*hStar_F_10d5m_v.*doesSvyIfAsked_I_10d5m_v.*wouldLieIfAsked_v);
PSVL_OO_NI_10d5m_5d1m_v = (1-r_v)*PSVL_NF_NI_10d5m_5d1m_v + r_v*mean((1 .-optsOutIfSees_OO_10d5m_v).*hStar_F_10d5m_v.*doesSvyIfAsked_NI_10d5m_v.*wouldLieIfAsked_5d1m_v);
PSVL_OO_I_10d5m_5d1m_v = (1-r_v)*PSVL_NF_I_10d5m_5d1m_v + r_v*mean((1 .-optsOutIfSees_OO_10d5m_v).*hStar_F_10d5m_v.*doesSvyIfAsked_I_10d5m_v.*wouldLieIfAsked_5d1m_v);

PSVL_OOV_NI_10d5m_v = (1-r_v)*PSVL_NF_NI_10d5m_v + r_v*mean((1 .-optsOutIfSees_OOV_10d5m_v).*hStar_FV_10d5m_v.*doesSvyIfAsked_I_10d5m_v.*wouldLieIfAsked_v);
PSVL_OOV_I_10d5m_v = (1-r_v)*PSVL_NF_I_10d5m_v + r_v*mean((1 .-optsOutIfSees_OOV_10d5m_v).*hStar_FV_10d5m_v.*doesSvyIfAsked_I_10d5m_v.*wouldLieIfAsked_v);
PSVL_OOV_NI_10d5m_5d1m_v = (1-r_v)*PSVL_NF_NI_10d5m_5d1m_v + r_v*mean((1 .-optsOutIfSees_OOV_10d5m_v).*hStar_FV_10d5m_v.*doesSvyIfAsked_I_10d5m_v.*wouldLieIfAsked_5d1m_v);
PSVL_OOV_I_10d5m_5d1m_v = (1-r_v)*PSVL_NF_I_10d5m_5d1m_v + r_v*mean((1 .-optsOutIfSees_OOV_10d5m_v).*hStar_FV_10d5m_v.*doesSvyIfAsked_I_10d5m_v.*wouldLieIfAsked_5d1m_v);

# !!!Non-Voters!!!

# PH = probability of being at home
PH_NF_0d5m_nv = h0_nv;
PH_F_0d5m_nv = (1-r_nv)*h0_nv + r_nv.*mean(hStar_F_0d5m_nv);
PH_FV_0d5m_nv = (1-r_nv)*h0_nv + r_nv.*mean(hStar_FV_0d5m_nv);
PH_OO_0d5m_nv = (1-r_nv)*h0_nv + r_nv.*mean((1 .-optsOutIfSees_OO_0d5m_nv).*hStar_F_0d5m_nv);
PH_OOV_0d5m_nv = (1-r_nv)*h0_nv + r_nv.*mean((1 .-optsOutIfSees_OOV_0d5m_nv).*hStar_FV_0d5m_nv);

PH_NF_10d10m_nv = h0_nv;
PH_F_10d10m_nv = (1-r_nv)*h0_nv + r_nv.*mean(hStar_F_10d10m_nv);
PH_FV_10d10m_nv = (1-r_nv)*h0_nv + r_nv.*mean(hStar_FV_10d10m_nv);
PH_OO_10d10m_nv = (1-r_nv)*h0_nv + r_nv.*mean((1 .-optsOutIfSees_OO_10d10m_nv).*hStar_F_10d10m_nv);
PH_OOV_10d10m_nv = (1-r_nv)*h0_nv + r_nv.*mean((1 .-optsOutIfSees_OOV_10d10m_nv).*hStar_FV_10d10m_nv);

PH_NF_10d5m_nv = h0_nv;
PH_F_10d5m_nv = (1-r_nv)*h0_nv + r_nv.*mean(hStar_F_10d5m_nv);
PH_FV_10d5m_nv = (1-r_nv)*h0_nv + r_nv.*mean(hStar_FV_10d5m_nv);
PH_OO_10d5m_nv = (1-r_nv)*h0_nv + r_nv.*mean((1 .-optsOutIfSees_OO_10d5m_nv).*hStar_F_10d5m_nv);
PH_OOV_10d5m_nv = (1-r_nv)*h0_nv + r_nv.*mean((1 .-optsOutIfSees_OOV_10d5m_nv).*hStar_FV_10d5m_nv);

# PSV=unconditional prob of doing the survey (not cond on opening door). PSV < PH mechanically

# 0d5m
PSV_NF_NI_0d5m_nv = h0_nv*mean(doesSvyIfAsked_NI_0d5m_nv);
PSV_NF_I_0d5m_nv = h0_nv*mean(doesSvyIfAsked_I_0d5m_nv);

PSV_F_NI_0d5m_nv = (1-r_nv)*PSV_NF_NI_0d5m_nv + r_nv*mean(hStar_F_0d5m_nv.*doesSvyIfAsked_NI_0d5m_nv);
PSV_F_I_0d5m_nv = (1-r_nv)*PSV_NF_I_0d5m_nv + r_nv*mean(hStar_F_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv);

PSV_FV_NI_0d5m_nv = (1-r_nv)*PSV_NF_NI_0d5m_nv + r_nv*mean(hStar_FV_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv);
PSV_FV_I_0d5m_nv = (1-r_nv)*PSV_NF_I_0d5m_nv + r_nv*mean(hStar_FV_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv);

PSV_OO_NI_0d5m_nv = (1-r_nv)*PSV_NF_NI_0d5m_nv + r_nv*mean((1 .-optsOutIfSees_OO_0d5m_nv).*hStar_F_0d5m_nv.*doesSvyIfAsked_NI_0d5m_nv);
PSV_OO_I_0d5m_nv = (1-r_nv)*PSV_NF_I_0d5m_nv + r_nv*mean((1 .-optsOutIfSees_OO_0d5m_nv).*hStar_F_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv);

PSV_OOV_NI_0d5m_nv = (1-r_nv)*PSV_NF_NI_0d5m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_0d5m_nv).*hStar_FV_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv);
PSV_OOV_I_0d5m_nv = (1-r_nv)*PSV_NF_I_0d5m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_0d5m_nv).*hStar_FV_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv);

# 10d10m
PSV_NF_NI_10d10m_nv = h0_nv*mean(doesSvyIfAsked_NI_10d10m_nv);
PSV_NF_I_10d10m_nv = h0_nv*mean(doesSvyIfAsked_I_10d10m_nv);

PSV_F_NI_10d10m_nv = (1-r_nv)*PSV_NF_NI_10d10m_nv + r_nv*mean(hStar_F_10d10m_nv.*doesSvyIfAsked_NI_10d10m_nv);
PSV_F_I_10d10m_nv = (1-r_nv)*PSV_NF_I_10d10m_nv + r_nv*mean(hStar_F_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv);

PSV_FV_NI_10d10m_nv = (1-r_nv)*PSV_NF_NI_10d10m_nv + r_nv*mean(hStar_FV_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv);
PSV_FV_I_10d10m_nv = (1-r_nv)*PSV_NF_I_10d10m_nv + r_nv*mean(hStar_FV_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv);

PSV_OO_NI_10d10m_nv = (1-r_nv)*PSV_NF_NI_10d10m_nv + r_nv*mean((1 .-optsOutIfSees_OO_10d10m_nv).*hStar_F_10d10m_nv.*doesSvyIfAsked_NI_10d10m_nv);
PSV_OO_I_10d10m_nv = (1-r_nv)*PSV_NF_I_10d10m_nv + r_nv*mean((1 .-optsOutIfSees_OO_10d10m_nv).*hStar_F_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv);

PSV_OOV_NI_10d10m_nv = (1-r_nv)*PSV_NF_NI_10d10m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_10d10m_nv).*hStar_FV_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv);
PSV_OOV_I_10d10m_nv = (1-r_nv)*PSV_NF_I_10d10m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_10d10m_nv).*hStar_FV_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv);

# 10d5m
PSV_NF_NI_10d5m_nv = h0_nv*mean(doesSvyIfAsked_NI_10d5m_nv);
PSV_NF_I_10d5m_nv = h0_nv*mean(doesSvyIfAsked_I_10d5m_nv);

PSV_F_NI_10d5m_nv = (1-r_nv)*PSV_NF_NI_10d5m_nv + r_nv*mean(hStar_F_10d5m_nv.*doesSvyIfAsked_NI_10d5m_nv);
PSV_F_I_10d5m_nv = (1-r_nv)*PSV_NF_I_10d5m_nv + r_nv*mean(hStar_F_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv);

PSV_FV_NI_10d5m_nv = (1-r_nv)*PSV_NF_NI_10d5m_nv + r_nv*mean(hStar_FV_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv);
PSV_FV_I_10d5m_nv = (1-r_nv)*PSV_NF_I_10d5m_nv + r_nv*mean(hStar_FV_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv);

PSV_OO_NI_10d5m_nv = (1-r_nv)*PSV_NF_NI_10d5m_nv + r_nv*mean((1 .-optsOutIfSees_OO_10d5m_nv).*hStar_F_10d5m_nv.*doesSvyIfAsked_NI_10d5m_nv);
PSV_OO_I_10d5m_nv = (1-r_nv)*PSV_NF_I_10d5m_nv + r_nv*mean((1 .-optsOutIfSees_OO_10d5m_nv).*hStar_F_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv);

PSV_OOV_NI_10d5m_nv = (1-r_nv)*PSV_NF_NI_10d5m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_10d5m_nv).*hStar_FV_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv);
PSV_OOV_I_10d5m_nv = (1-r_nv)*PSV_NF_I_10d5m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_10d5m_nv).*hStar_FV_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv);

# POO=prob of opting out (not conditional on seeing flyer). Scaled by baseline likelihood of being at home
POO_OO_0d5m_nv =  h0_nv*r_nv*mean(optsOutIfSees_OO_0d5m_nv);
POO_OOV_0d5m_nv = h0_nv*r_nv*mean(optsOutIfSees_OOV_0d5m_nv);

POO_OO_10d10m_nv =  h0_nv*r_nv*mean(optsOutIfSees_OO_10d10m_nv);
POO_OOV_10d10m_nv = h0_nv*r_nv*mean(optsOutIfSees_OOV_10d10m_nv);

POO_OO_10d5m_nv =  h0_nv*r_nv*mean(optsOutIfSees_OO_10d5m_nv);
POO_OOV_10d5m_nv = h0_nv*r_nv*mean(optsOutIfSees_OOV_10d5m_nv);

# PSVL = unconditional percent who do survey and lie No flyer treatment only, simplifies later code

# PL=cond on agreeing to do the survey, did you lie? Incentive to lie is a surprise later (doesn't affect PH or PSV)

# 0d5m, 5d1m incentive
PSVL_NF_NI_0d5m_nv = mean(h0_nv.*doesSvyIfAsked_NI_0d5m_nv.*wouldLieIfAsked_nv);
PSVL_NF_I_0d5m_nv = mean(h0_nv.*doesSvyIfAsked_I_0d5m_nv.*wouldLieIfAsked_nv);
PSVL_NF_NI_0d5m_5d1m_nv = mean(h0_nv.*doesSvyIfAsked_NI_0d5m_nv.*wouldLieIfAsked_5d1m_nv);
PSVL_NF_I_0d5m_5d1m_nv = mean(h0_nv.*doesSvyIfAsked_I_0d5m_nv.*wouldLieIfAsked_5d1m_nv);

PSVL_F_NI_0d5m_nv = (1-r_nv)*PSVL_NF_NI_0d5m_nv + r_nv*mean(hStar_F_0d5m_nv.*doesSvyIfAsked_NI_0d5m_nv.*wouldLieIfAsked_nv);
PSVL_F_I_0d5m_nv = (1-r_nv)*PSVL_NF_I_0d5m_nv + r_nv*mean(hStar_F_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv.*wouldLieIfAsked_nv);
PSVL_F_NI_0d5m_5d1m_nv = (1-r_nv)*PSVL_NF_NI_0d5m_5d1m_nv + r_nv*mean(hStar_F_0d5m_nv.*doesSvyIfAsked_NI_0d5m_nv.*wouldLieIfAsked_5d1m_nv);
PSVL_F_I_0d5m_5d1m_nv = (1-r_nv)*PSVL_NF_I_0d5m_5d1m_nv + r_nv*mean(hStar_F_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv.*wouldLieIfAsked_5d1m_nv);

PSVL_FV_NI_0d5m_nv = (1-r_nv)*PSVL_NF_NI_0d5m_nv + r_nv*mean(hStar_FV_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv.*wouldLieIfAsked_nv);
PSVL_FV_I_0d5m_nv = (1-r_nv)*PSVL_NF_I_0d5m_nv + r_nv*mean(hStar_FV_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv.*wouldLieIfAsked_nv);
PSVL_FV_NI_0d5m_5d1m_nv = (1-r_nv)*PSVL_NF_NI_0d5m_5d1m_nv + r_nv*mean(hStar_FV_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv.*wouldLieIfAsked_5d1m_nv);
PSVL_FV_I_0d5m_5d1m_nv = (1-r_nv)*PSVL_NF_I_0d5m_5d1m_nv + r_nv*mean(hStar_FV_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv.*wouldLieIfAsked_5d1m_nv);

PSVL_OO_NI_0d5m_nv = (1-r_nv)*PSVL_NF_NI_0d5m_nv + r_nv*mean((1 .-optsOutIfSees_OO_0d5m_nv).*hStar_F_0d5m_nv.*doesSvyIfAsked_NI_0d5m_nv.*wouldLieIfAsked_nv);
PSVL_OO_I_0d5m_nv = (1-r_nv)*PSVL_NF_I_0d5m_nv + r_nv*mean((1 .-optsOutIfSees_OO_0d5m_nv).*hStar_F_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv.*wouldLieIfAsked_nv);
PSVL_OO_NI_0d5m_5d1m_nv = (1-r_nv)*PSVL_NF_NI_0d5m_5d1m_nv + r_nv*mean((1 .-optsOutIfSees_OO_0d5m_nv).*hStar_F_0d5m_nv.*doesSvyIfAsked_NI_0d5m_nv.*wouldLieIfAsked_5d1m_nv);
PSVL_OO_I_0d5m_5d1m_nv = (1-r_nv)*PSVL_NF_I_0d5m_5d1m_nv + r_nv*mean((1 .-optsOutIfSees_OO_0d5m_nv).*hStar_F_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv.*wouldLieIfAsked_5d1m_nv);

PSVL_OOV_NI_0d5m_nv = (1-r_nv)*PSVL_NF_NI_0d5m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_0d5m_nv).*hStar_FV_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv.*wouldLieIfAsked_nv);
PSVL_OOV_I_0d5m_nv = (1-r_nv)*PSVL_NF_I_0d5m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_0d5m_nv).*hStar_FV_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv.*wouldLieIfAsked_nv);
PSVL_OOV_NI_0d5m_5d1m_nv = (1-r_nv)*PSVL_NF_NI_0d5m_5d1m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_0d5m_nv).*hStar_FV_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv.*wouldLieIfAsked_5d1m_nv);
PSVL_OOV_I_0d5m_5d1m_nv = (1-r_nv)*PSVL_NF_I_0d5m_5d1m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_0d5m_nv).*hStar_FV_0d5m_nv.*doesSvyIfAsked_I_0d5m_nv.*wouldLieIfAsked_5d1m_nv);

# 10d10m, 8m incentive
PSVL_NF_NI_10d10m_nv = mean(h0_nv.*doesSvyIfAsked_NI_10d10m_nv.*wouldLieIfAsked_nv);
PSVL_NF_I_10d10m_nv = mean(h0_nv.*doesSvyIfAsked_I_10d10m_nv.*wouldLieIfAsked_nv);
PSVL_NF_NI_10d10m_8m_nv = mean(h0_nv.*doesSvyIfAsked_NI_10d10m_nv.*wouldLieIfAsked_8m_nv);
PSVL_NF_I_10d10m_8m_nv = mean(h0_nv.*doesSvyIfAsked_I_10d10m_nv.*wouldLieIfAsked_8m_nv);

PSVL_F_NI_10d10m_nv = (1-r_nv)*PSVL_NF_NI_10d10m_nv + r_nv*mean(hStar_F_10d10m_nv.*doesSvyIfAsked_NI_10d10m_nv.*wouldLieIfAsked_nv);
PSVL_F_I_10d10m_nv = (1-r_nv)*PSVL_NF_I_10d10m_nv + r_nv*mean(hStar_F_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv.*wouldLieIfAsked_nv);
PSVL_F_NI_10d10m_8m_nv = (1-r_nv)*PSVL_NF_NI_10d10m_8m_nv + r_nv*mean(hStar_F_10d10m_nv.*doesSvyIfAsked_NI_10d10m_nv.*wouldLieIfAsked_8m_nv);
PSVL_F_I_10d10m_8m_nv = (1-r_nv)*PSVL_NF_I_10d10m_8m_nv + r_nv*mean(hStar_F_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv.*wouldLieIfAsked_8m_nv);

PSVL_FV_NI_10d10m_nv = (1-r_nv)*PSVL_NF_NI_10d10m_nv + r_nv*mean(hStar_FV_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv.*wouldLieIfAsked_nv);
PSVL_FV_I_10d10m_nv = (1-r_nv)*PSVL_NF_I_10d10m_nv + r_nv*mean(hStar_FV_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv.*wouldLieIfAsked_nv);
PSVL_FV_NI_10d10m_8m_nv = (1-r_nv)*PSVL_NF_NI_10d10m_8m_nv + r_nv*mean(hStar_FV_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv.*wouldLieIfAsked_8m_nv);
PSVL_FV_I_10d10m_8m_nv = (1-r_nv)*PSVL_NF_I_10d10m_8m_nv + r_nv*mean(hStar_FV_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv.*wouldLieIfAsked_8m_nv);

PSVL_OO_NI_10d10m_nv = (1-r_nv)*PSVL_NF_NI_10d10m_nv + r_nv*mean((1 .-optsOutIfSees_OO_10d10m_nv).*hStar_F_10d10m_nv.*doesSvyIfAsked_NI_10d10m_nv.*wouldLieIfAsked_nv);
PSVL_OO_I_10d10m_nv = (1-r_nv)*PSVL_NF_I_10d10m_nv + r_nv*mean((1 .-optsOutIfSees_OO_10d10m_nv).*hStar_F_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv.*wouldLieIfAsked_nv);
PSVL_OO_NI_10d10m_8m_nv = (1-r_nv)*PSVL_NF_NI_10d10m_8m_nv + r_nv*mean((1 .-optsOutIfSees_OO_10d10m_nv).*hStar_F_10d10m_nv.*doesSvyIfAsked_NI_10d10m_nv.*wouldLieIfAsked_8m_nv);
PSVL_OO_I_10d10m_8m_nv = (1-r_nv)*PSVL_NF_I_10d10m_8m_nv + r_nv*mean((1 .-optsOutIfSees_OO_10d10m_nv).*hStar_F_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv.*wouldLieIfAsked_8m_nv);

PSVL_OOV_NI_10d10m_nv = (1-r_nv)*PSVL_NF_NI_10d10m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_10d10m_nv).*hStar_FV_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv.*wouldLieIfAsked_nv);
PSVL_OOV_I_10d10m_nv = (1-r_nv)*PSVL_NF_I_10d10m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_10d10m_nv).*hStar_FV_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv.*wouldLieIfAsked_nv);
PSVL_OOV_NI_10d10m_8m_nv = (1-r_nv)*PSVL_NF_NI_10d10m_8m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_10d10m_nv).*hStar_FV_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv.*wouldLieIfAsked_8m_nv);
PSVL_OOV_I_10d10m_8m_nv = (1-r_nv)*PSVL_NF_I_10d10m_8m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_10d10m_nv).*hStar_FV_10d10m_nv.*doesSvyIfAsked_I_10d10m_nv.*wouldLieIfAsked_8m_nv);

# 10d5m, 5d1m incentive
PSVL_NF_NI_10d5m_nv = mean(h0_nv.*doesSvyIfAsked_NI_10d5m_nv.*wouldLieIfAsked_nv);
PSVL_NF_I_10d5m_nv = mean(h0_nv.*doesSvyIfAsked_I_10d5m_nv.*wouldLieIfAsked_nv);
PSVL_NF_NI_10d5m_5d1m_nv = mean(h0_nv.*doesSvyIfAsked_NI_10d5m_nv.*wouldLieIfAsked_5d1m_nv);
PSVL_NF_I_10d5m_5d1m_nv = mean(h0_nv.*doesSvyIfAsked_I_10d5m_nv.*wouldLieIfAsked_5d1m_nv);

PSVL_F_NI_10d5m_nv = (1-r_nv)*PSVL_NF_NI_10d5m_nv + r_nv*mean(hStar_F_10d5m_nv.*doesSvyIfAsked_NI_10d5m_nv.*wouldLieIfAsked_nv);
PSVL_F_I_10d5m_nv = (1-r_nv)*PSVL_NF_I_10d5m_nv + r_nv*mean(hStar_F_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv.*wouldLieIfAsked_nv);
PSVL_F_NI_10d5m_5d1m_nv = (1-r_nv)*PSVL_NF_NI_10d5m_5d1m_nv + r_nv*mean(hStar_F_10d5m_nv.*doesSvyIfAsked_NI_10d5m_nv.*wouldLieIfAsked_5d1m_nv);
PSVL_F_I_10d5m_5d1m_nv = (1-r_nv)*PSVL_NF_I_10d5m_5d1m_nv + r_nv*mean(hStar_F_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv.*wouldLieIfAsked_5d1m_nv);

PSVL_FV_NI_10d5m_nv = (1-r_nv)*PSVL_NF_NI_10d5m_nv + r_nv*mean(hStar_FV_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv.*wouldLieIfAsked_nv);
PSVL_FV_I_10d5m_nv = (1-r_nv)*PSVL_NF_I_10d5m_nv + r_nv*mean(hStar_FV_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv.*wouldLieIfAsked_nv);
PSVL_FV_NI_10d5m_5d1m_nv = (1-r_nv)*PSVL_NF_NI_10d5m_5d1m_nv + r_nv*mean(hStar_FV_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv.*wouldLieIfAsked_5d1m_nv);
PSVL_FV_I_10d5m_5d1m_nv = (1-r_nv)*PSVL_NF_I_10d5m_5d1m_nv + r_nv*mean(hStar_FV_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv.*wouldLieIfAsked_5d1m_nv);

PSVL_OO_NI_10d5m_nv = (1-r_nv)*PSVL_NF_NI_10d5m_nv + r_nv*mean((1 .-optsOutIfSees_OO_10d5m_nv).*hStar_F_10d5m_nv.*doesSvyIfAsked_NI_10d5m_nv.*wouldLieIfAsked_nv);
PSVL_OO_I_10d5m_nv = (1-r_nv)*PSVL_NF_I_10d5m_nv + r_nv*mean((1 .-optsOutIfSees_OO_10d5m_nv).*hStar_F_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv.*wouldLieIfAsked_nv);
PSVL_OO_NI_10d5m_5d1m_nv = (1-r_nv)*PSVL_NF_NI_10d5m_5d1m_nv + r_nv*mean((1 .-optsOutIfSees_OO_10d5m_nv).*hStar_F_10d5m_nv.*doesSvyIfAsked_NI_10d5m_nv.*wouldLieIfAsked_5d1m_nv);
PSVL_OO_I_10d5m_5d1m_nv = (1-r_nv)*PSVL_NF_I_10d5m_5d1m_nv + r_nv*mean((1 .-optsOutIfSees_OO_10d5m_nv).*hStar_F_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv.*wouldLieIfAsked_5d1m_nv);

PSVL_OOV_NI_10d5m_nv = (1-r_nv)*PSVL_NF_NI_10d5m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_10d5m_nv).*hStar_FV_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv.*wouldLieIfAsked_nv);
PSVL_OOV_I_10d5m_nv = (1-r_nv)*PSVL_NF_I_10d5m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_10d5m_nv).*hStar_FV_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv.*wouldLieIfAsked_nv);
PSVL_OOV_NI_10d5m_5d1m_nv = (1-r_nv)*PSVL_NF_NI_10d5m_5d1m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_10d5m_nv).*hStar_FV_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv.*wouldLieIfAsked_5d1m_nv);
PSVL_OOV_I_10d5m_5d1m_nv = (1-r_nv)*PSVL_NF_I_10d5m_5d1m_nv + r_nv*mean((1 .-optsOutIfSees_OOV_10d5m_nv).*hStar_FV_10d5m_nv.*doesSvyIfAsked_I_10d5m_nv.*wouldLieIfAsked_5d1m_nv);

# 30 PH moments (5 tx * 3 length * 2 v/nv)
    sm = Vector{Real}(undef, 101)
    sm[1]	=		PH_NF_0d5m_v
    sm[2]	=		PH_NF_10d10m_v
    sm[3]	=		PH_NF_10d5m_v
    sm[4]	=		PH_F_0d5m_v
    sm[5]	=		PH_F_10d10m_v
    sm[6]	=		PH_F_10d5m_v
    sm[7]	=		PH_FV_0d5m_v
    sm[8]	=		PH_FV_10d10m_v
    sm[9]	=		PH_FV_10d5m_v
    sm[10]	=		PH_OO_0d5m_v
    sm[11]	=		PH_OO_10d10m_v
    sm[12]	=		PH_OO_10d5m_v
    sm[13]	=		PH_OOV_0d5m_v
    sm[14]	=		PH_OOV_10d10m_v
    sm[15]	=		PH_OOV_10d5m_v
    sm[16]	=		PH_NF_0d5m_nv
    sm[17]	=		PH_NF_10d10m_nv
    sm[18]	=		PH_NF_10d5m_nv
    sm[19]	=		PH_F_0d5m_nv
    sm[20]	=		PH_F_10d10m_nv
    sm[21]	=		PH_F_10d5m_nv
    sm[22]	=		PH_FV_0d5m_nv
    sm[23]	=		PH_FV_10d10m_nv
    sm[24]	=		PH_FV_10d5m_nv
    sm[25]	=		PH_OO_0d5m_nv
    sm[26]	=		PH_OO_10d10m_nv
    sm[27]	=		PH_OO_10d5m_nv
    sm[28]	=		PH_OOV_0d5m_nv
    sm[29]	=		PH_OOV_10d10m_nv
    sm[30]	=		PH_OOV_10d5m_nv

# 30 PSV moments (5 tx * 3 length * 2 v/nv). Taking 50/50 average across NI and I treatments
    sm[31]  =   mean([  PSV_NF_I_0d5m_v,     PSV_NF_NI_0d5m_v    ])
    sm[32]  =   mean([  PSV_NF_I_10d10m_v,   PSV_NF_NI_10d10m_v  ])
    sm[33]  =   mean([  PSV_NF_I_10d5m_v,    PSV_NF_NI_10d5m_v   ])
    sm[34]  =   mean([  PSV_F_I_0d5m_v,      PSV_F_NI_0d5m_v     ])
    sm[35]  =   mean([  PSV_F_I_10d10m_v,    PSV_F_NI_10d10m_v   ])
    sm[36]  =   mean([  PSV_F_I_10d5m_v,     PSV_F_NI_10d5m_v    ])
    sm[37]  =   mean([  PSV_FV_I_0d5m_v,     PSV_FV_NI_0d5m_v    ])
    sm[38]  =   mean([  PSV_FV_I_10d10m_v,   PSV_FV_NI_10d10m_v  ])
    sm[39]  =   mean([  PSV_FV_I_10d5m_v,    PSV_FV_NI_10d5m_v   ])
    sm[40] =   mean([  PSV_OO_I_0d5m_v,     PSV_OO_NI_0d5m_v    ])
    sm[41] =   mean([  PSV_OO_I_10d10m_v,   PSV_OO_NI_10d10m_v  ])
    sm[42] =   mean([  PSV_OO_I_10d5m_v,    PSV_OO_NI_10d5m_v   ])
    sm[43] =   mean([  PSV_OOV_I_0d5m_v,    PSV_OOV_NI_0d5m_v   ])
    sm[44] =   mean([  PSV_OOV_I_10d10m_v,  PSV_OOV_NI_10d10m_v ])
    sm[45] =   mean([  PSV_OOV_I_10d5m_v,   PSV_OOV_NI_10d5m_v  ])
    sm[46] =   mean([  PSV_NF_I_0d5m_nv,    PSV_NF_NI_0d5m_nv   ])
    sm[47] =   mean([  PSV_NF_I_10d10m_nv,  PSV_NF_NI_10d10m_nv ])
    sm[48] =   mean([  PSV_NF_I_10d5m_nv,   PSV_NF_NI_10d5m_nv  ])
    sm[49] =   mean([  PSV_F_I_0d5m_nv,     PSV_F_NI_0d5m_nv    ])
    sm[50] =   mean([  PSV_F_I_10d10m_nv,   PSV_F_NI_10d10m_nv  ])
    sm[51] =   mean([  PSV_F_I_10d5m_nv,    PSV_F_NI_10d5m_nv   ])
    sm[52] =   mean([  PSV_FV_I_0d5m_nv,    PSV_FV_NI_0d5m_nv   ])
    sm[53] =   mean([  PSV_FV_I_10d10m_nv,  PSV_FV_NI_10d10m_nv ])
    sm[54] =   mean([  PSV_FV_I_10d5m_nv,   PSV_FV_NI_10d5m_nv  ])
    sm[55] =   mean([  PSV_OO_I_0d5m_nv,    PSV_OO_NI_0d5m_nv   ])
    sm[56] =   mean([  PSV_OO_I_10d10m_nv,  PSV_OO_NI_10d10m_nv ])
    sm[57] =   mean([  PSV_OO_I_10d5m_nv,   PSV_OO_NI_10d5m_nv  ])
    sm[58] =   mean([  PSV_OOV_I_0d5m_nv,   PSV_OOV_NI_0d5m_nv  ])
    sm[59] =   mean([  PSV_OOV_I_10d10m_nv, PSV_OOV_NI_10d10m_nv])
    sm[60] =   mean([  PSV_OOV_I_10d5m_nv,  PSV_OOV_NI_10d5m_nv ])

# 12 POO moments (2 tx * 3 length * 2 v/nv)
    sm[61]	=		POO_OO_0d5m_v
    sm[62]	=		POO_OO_10d10m_v
    sm[63]	=		POO_OO_10d5m_v
    sm[64]	=		POO_OOV_0d5m_v
    sm[65]	=		POO_OOV_10d10m_v
    sm[66]	=		POO_OOV_10d5m_v
    sm[67]	=		POO_OO_0d5m_nv
    sm[68]	=		POO_OO_10d10m_nv
    sm[69]	=		POO_OO_10d5m_nv
    sm[70]	=		POO_OOV_0d5m_nv
    sm[71]	=		POO_OOV_10d10m_nv
    sm[72]	=		POO_OOV_10d5m_nv

# 20 PSV by info moments (5 tx * 2 I/NI * 2 v/nv)
    sm[73]   =   mean([  PSV_NF_NI_0d5m_v,    PSV_NF_NI_10d10m_v,      PSV_NF_NI_10d5m_v   ])
    sm[74]   =   mean([  PSV_NF_I_0d5m_v,     PSV_NF_I_10d10m_v,       PSV_NF_I_10d5m_v  ])
    sm[75]   =   mean([  PSV_F_NI_0d5m_v,     PSV_F_NI_10d10m_v,       PSV_F_NI_10d5m_v    ])
    sm[76]   =   mean([  PSV_F_I_0d5m_v,      PSV_F_I_10d10m_v,       PSV_F_I_10d5m_v   ])
    sm[77]   =   mean([  PSV_FV_NI_0d5m_v,    PSV_FV_NI_10d10m_v,      PSV_FV_NI_10d5m_v   ])
    sm[78]   =   mean([  PSV_FV_I_0d5m_v,     PSV_FV_I_10d10m_v,       PSV_FV_I_10d5m_v  ])
    sm[79]   =   mean([  PSV_OO_NI_0d5m_v,    PSV_OO_NI_10d10m_v,      PSV_OO_NI_10d5m_v   ])
    sm[80]   =   mean([  PSV_OO_I_0d5m_v,     PSV_OO_I_10d10m_v,       PSV_OO_I_10d5m_v  ])
    sm[81]   =   mean([  PSV_OOV_NI_0d5m_v,   PSV_OOV_NI_10d10m_v,     PSV_OOV_NI_10d5m_v  ])
    sm[82]   =  mean([  PSV_OOV_I_0d5m_v,    PSV_OOV_I_10d10m_v,      PSV_OOV_I_10d5m_v ])
    sm[83]  =   mean([  PSV_NF_NI_0d5m_nv,   PSV_NF_NI_10d10m_nv,     PSV_NF_NI_10d5m_nv  ])
    sm[84]  =   mean([  PSV_NF_I_0d5m_nv,    PSV_NF_I_10d10m_nv,      PSV_NF_I_10d5m_nv ])
    sm[85]  =   mean([  PSV_F_NI_0d5m_nv,    PSV_F_NI_10d10m_nv,      PSV_F_NI_10d5m_nv   ])
    sm[86]  =   mean([  PSV_F_I_0d5m_nv,     PSV_F_I_10d10m_nv,       PSV_F_I_10d5m_nv  ])
    sm[87]  =   mean([  PSV_FV_NI_0d5m_nv,   PSV_FV_NI_10d10m_nv,     PSV_FV_NI_10d5m_nv  ])
    sm[88]  =   mean([  PSV_FV_I_0d5m_nv,    PSV_FV_I_10d10m_nv,      PSV_FV_I_10d5m_nv ])
    sm[89]  =   mean([  PSV_OO_NI_0d5m_nv,   PSV_OO_NI_10d10m_nv,     PSV_OO_NI_10d5m_nv  ])
    sm[90]  =   mean([  PSV_OO_I_0d5m_nv,    PSV_OO_I_10d10m_nv,      PSV_OO_I_10d5m_nv ])
    sm[91]  =   mean([  PSV_OOV_NI_0d5m_nv,  PSV_OOV_NI_10d10m_nv,    PSV_OOV_NI_10d5m_nv ])
    sm[92]  =   mean([  PSV_OOV_I_0d5m_nv,   PSV_OOV_I_10d10m_nv,     PSV_OOV_I_10d5m_nv ])

# 8 PL moments (1 tx * 2 10m/5m * 2 incentives)

# Empirical moments are sum of people lying in relevant tx divided by the sum of people answering the survey in relevant tx.
    sm[93] = mean([PSVL_NF_NI_0d5m_v, PSVL_NF_I_0d5m_v, PSVL_NF_NI_10d5m_v, PSVL_NF_I_10d5m_v,
                    PSVL_F_NI_0d5m_v, PSVL_F_I_0d5m_v, PSVL_F_NI_10d5m_v, PSVL_F_I_10d5m_v,
                    PSVL_FV_NI_0d5m_v, PSVL_FV_I_0d5m_v, PSVL_FV_NI_10d5m_v, PSVL_FV_I_10d5m_v,
                    PSVL_OO_NI_0d5m_v, PSVL_OO_I_0d5m_v, PSVL_OO_NI_10d5m_v, PSVL_OO_I_10d5m_v,
                    PSVL_OOV_NI_0d5m_v, PSVL_OOV_I_0d5m_v, PSVL_OOV_NI_10d5m_v, PSVL_OOV_I_10d5m_v])/mean([PSV_NF_NI_0d5m_v,
                    PSV_NF_I_0d5m_v, PSV_NF_NI_10d5m_v, PSV_NF_I_10d5m_v,
                    PSV_F_NI_0d5m_v, PSV_F_I_0d5m_v, PSV_F_NI_10d5m_v, PSV_F_I_10d5m_v,
                    PSV_FV_NI_0d5m_v, PSV_FV_I_0d5m_v, PSV_FV_NI_10d5m_v, PSV_FV_I_10d5m_v,
                    PSV_OO_NI_0d5m_v, PSV_OO_I_0d5m_v,  PSV_OO_NI_10d5m_v, PSV_OO_I_10d5m_v,
                    PSV_OOV_NI_0d5m_v, PSV_OOV_I_0d5m_v, PSV_OOV_NI_10d5m_v, PSV_OOV_I_10d5m_v])

    sm[94] = mean([PSVL_NF_NI_0d5m_5d1m_v, PSVL_NF_I_0d5m_5d1m_v, PSVL_NF_NI_10d5m_5d1m_v, PSVL_NF_I_10d5m_5d1m_v,
                    PSVL_F_NI_0d5m_5d1m_v, PSVL_F_I_0d5m_5d1m_v, PSVL_F_NI_10d5m_5d1m_v, PSVL_F_I_10d5m_5d1m_v,
                    PSVL_FV_NI_0d5m_5d1m_v, PSVL_FV_I_0d5m_5d1m_v, PSVL_FV_NI_10d5m_5d1m_v, PSVL_FV_I_10d5m_5d1m_v,
                    PSVL_OO_NI_0d5m_5d1m_v, PSVL_OO_I_0d5m_5d1m_v, PSVL_OO_NI_10d5m_5d1m_v, PSVL_OO_I_10d5m_5d1m_v,
                    PSVL_OOV_NI_0d5m_5d1m_v, PSVL_OOV_I_0d5m_5d1m_v, PSVL_OOV_NI_10d5m_5d1m_v, PSVL_OOV_I_10d5m_5d1m_v])/ mean([PSV_NF_NI_0d5m_v,
                    PSV_NF_I_0d5m_v, PSV_NF_NI_10d5m_v, PSV_NF_I_10d5m_v,
                    PSV_F_NI_0d5m_v, PSV_F_I_0d5m_v, PSV_F_NI_10d5m_v, PSV_F_I_10d5m_v,
                    PSV_FV_NI_0d5m_v, PSV_FV_I_0d5m_v, PSV_FV_NI_10d5m_v, PSV_FV_I_10d5m_v,
                    PSV_OO_NI_0d5m_v, PSV_OO_I_0d5m_v, PSV_OO_NI_10d5m_v, PSV_OO_I_10d5m_v,
                    PSV_OOV_NI_0d5m_v, PSV_OOV_I_0d5m_v, PSV_OOV_NI_10d5m_v, PSV_OOV_I_10d5m_v])

    sm[95] = mean([PSVL_NF_NI_10d10m_v, PSVL_NF_I_10d10m_v,
                    PSVL_F_NI_10d10m_v, PSVL_F_I_10d10m_v,
                    PSVL_FV_NI_10d10m_v, PSVL_FV_I_10d10m_v,
                    PSVL_OO_NI_10d10m_v, PSVL_OO_I_10d10m_v,
                    PSVL_OOV_NI_10d10m_v, PSVL_OOV_I_10d10m_v])/ mean([PSV_NF_NI_10d10m_v,
                    PSV_NF_I_10d10m_v,
                    PSV_F_NI_10d10m_v, PSV_F_I_10d10m_v,
                    PSV_FV_NI_10d10m_v, PSV_FV_I_10d10m_v,
                    PSV_OO_NI_10d10m_v, PSV_OO_I_10d10m_v,
                    PSV_OOV_NI_10d10m_v, PSV_OOV_I_10d10m_v])

    sm[96] = mean([PSVL_NF_NI_10d10m_8m_v, PSVL_NF_I_10d10m_8m_v,
                    PSVL_F_NI_10d10m_8m_v, PSVL_F_I_10d10m_8m_v,
                    PSVL_FV_NI_10d10m_8m_v, PSVL_FV_I_10d10m_8m_v,
                    PSVL_OO_NI_10d10m_8m_v, PSVL_OO_I_10d10m_8m_v,
                    PSVL_OOV_NI_10d10m_8m_v, PSVL_OOV_I_10d10m_8m_v])/ mean([PSV_NF_NI_10d10m_v,
                    PSV_NF_I_10d10m_v,
                    PSV_F_NI_10d10m_v, PSV_F_I_10d10m_v, PSV_FV_NI_10d10m_v, PSV_FV_I_10d10m_v,
                    PSV_OO_NI_10d10m_v, PSV_OO_I_10d10m_v,
                    PSV_OOV_NI_10d10m_v, PSV_OOV_I_10d10m_v])

    sm[97] = mean([PSVL_NF_NI_0d5m_nv, PSVL_NF_I_0d5m_nv, PSVL_NF_NI_10d5m_nv, PSVL_NF_I_10d5m_nv,
                    PSVL_F_NI_0d5m_nv, PSVL_F_I_0d5m_nv, PSVL_F_NI_10d5m_nv, PSVL_F_I_10d5m_nv,
                    PSVL_FV_NI_0d5m_nv, PSVL_FV_I_0d5m_nv, PSVL_FV_NI_10d5m_nv, PSVL_FV_I_10d5m_nv,
                    PSVL_OO_NI_0d5m_nv, PSVL_OO_I_0d5m_nv, PSVL_OO_NI_10d5m_nv, PSVL_OO_I_10d5m_nv,
                    PSVL_OOV_NI_0d5m_nv, PSVL_OOV_I_0d5m_nv, PSVL_OOV_NI_10d5m_nv, PSVL_OOV_I_10d5m_nv])/ mean([PSV_NF_NI_0d5m_nv,
                    PSV_NF_I_0d5m_nv, PSV_NF_NI_10d5m_nv, PSV_NF_I_10d5m_nv,
                    PSV_F_NI_0d5m_nv, PSV_F_I_0d5m_nv, PSV_F_NI_10d5m_nv, PSV_F_I_10d5m_nv,
                    PSV_FV_NI_0d5m_nv, PSV_FV_I_0d5m_nv, PSV_FV_NI_10d5m_nv, PSV_FV_I_10d5m_nv,
                    PSV_OO_NI_0d5m_nv, PSV_OO_I_0d5m_nv,  PSV_OO_NI_10d5m_nv, PSV_OO_I_10d5m_nv,
                    PSV_OOV_NI_0d5m_nv, PSV_OOV_I_0d5m_nv, PSV_OOV_NI_10d5m_nv, PSV_OOV_I_10d5m_nv])

    sm[98] = mean([PSVL_NF_NI_0d5m_5d1m_nv, PSVL_NF_I_0d5m_5d1m_nv, PSVL_NF_NI_10d5m_5d1m_nv, PSVL_NF_I_10d5m_5d1m_nv, PSVL_F_NI_0d5m_5d1m_nv, PSVL_F_I_0d5m_5d1m_nv, PSVL_F_NI_10d5m_5d1m_nv, PSVL_F_I_10d5m_5d1m_nv, PSVL_FV_NI_0d5m_5d1m_nv, PSVL_FV_I_0d5m_5d1m_nv,PSVL_FV_NI_10d5m_5d1m_nv, PSVL_FV_I_10d5m_5d1m_nv,PSVL_OO_NI_0d5m_5d1m_nv, PSVL_OO_I_0d5m_5d1m_nv, PSVL_OO_NI_10d5m_5d1m_nv, PSVL_OO_I_10d5m_5d1m_nv,PSVL_OOV_NI_0d5m_5d1m_nv, PSVL_OOV_I_0d5m_5d1m_nv, PSVL_OOV_NI_10d5m_5d1m_nv, PSVL_OOV_I_10d5m_5d1m_nv ])/mean([PSV_NF_NI_0d5m_nv, PSV_NF_I_0d5m_nv, PSV_NF_NI_10d5m_nv, PSV_NF_I_10d5m_nv,
PSV_F_NI_0d5m_nv, PSV_F_I_0d5m_nv, PSV_F_NI_10d5m_nv, PSV_F_I_10d5m_nv,
PSV_FV_NI_0d5m_nv, PSV_FV_I_0d5m_nv, PSV_FV_NI_10d5m_nv, PSV_FV_I_10d5m_nv,
PSV_OO_NI_0d5m_nv, PSV_OO_I_0d5m_nv,  PSV_OO_NI_10d5m_nv, PSV_OO_I_10d5m_nv,
PSV_OOV_NI_0d5m_nv, PSV_OOV_I_0d5m_nv, PSV_OOV_NI_10d5m_nv, PSV_OOV_I_10d5m_nv])

    sm[99] =	mean([PSVL_NF_NI_10d10m_nv, PSVL_NF_I_10d10m_nv,
                    PSVL_F_NI_10d10m_nv, PSVL_F_I_10d10m_nv,
                    PSVL_FV_NI_10d10m_nv, PSVL_FV_I_10d10m_nv,
                    PSVL_OO_NI_10d10m_nv, PSVL_OO_I_10d10m_nv,
                    PSVL_OOV_NI_10d10m_nv, PSVL_OOV_I_10d10m_nv])/mean([PSV_NF_NI_10d10m_nv,
                    PSV_NF_I_10d10m_nv,
                    PSV_F_NI_10d10m_nv, PSV_F_I_10d10m_nv,
                    PSV_FV_NI_10d10m_nv, PSV_FV_I_10d10m_nv,
                    PSV_OO_NI_10d10m_nv, PSV_OO_I_10d10m_nv,
                    PSV_OOV_NI_10d10m_nv, PSV_OOV_I_10d10m_nv])

    sm[100] =	mean([PSVL_NF_NI_10d10m_8m_nv, PSVL_NF_I_10d10m_8m_nv,
                    PSVL_F_NI_10d10m_8m_nv, PSVL_F_I_10d10m_8m_nv,
                    PSVL_FV_NI_10d10m_8m_nv, PSVL_FV_I_10d10m_8m_nv,
                    PSVL_OO_NI_10d10m_8m_nv, PSVL_OO_I_10d10m_8m_nv,
                    PSVL_OOV_NI_10d10m_8m_nv, PSVL_OOV_I_10d10m_8m_nv])/mean([PSV_NF_NI_10d10m_nv,
                    PSV_NF_I_10d10m_nv,
                    PSV_F_NI_10d10m_nv, PSV_F_I_10d10m_nv,
                    PSV_FV_NI_10d10m_nv, PSV_FV_I_10d10m_nv,
                    PSV_OO_NI_10d10m_nv, PSV_OO_I_10d10m_nv,
                    PSV_OOV_NI_10d10m_nv, PSV_OOV_I_10d10m_nv])

	sm[101] = Turnout_control     # baseline turnout

    return sm
end

# ╔═╡ c6468ee0-3588-11ec-0d41-cbfce972d21c
begin

# 3. Estimation

# Define the number N of simulated individuals and the draws from a standard normal distribution to compute the simulated s, s_v, s_nv, eps in the function

sim_voters = 750000
rand1 = rand(Normal(0.0, 1.0), sim_voters)
rand2 = rand(Normal(0.0, 1.0), sim_voters)
rand3 = rand(Normal(0.0, 1.0), sim_voters)
rand4 = rand(Normal(0.0, 1.0), sim_voters)
rand_vec = [rand1,rand2,rand3,rand4]

end;

# ╔═╡ c616cc50-3588-11ec-25bc-1b52b9c8694d
# This is the function we want to minimize: the weighted sum of the squared differences between empirical and simulated moments

function criterion(parameters, rand_set)

    simMoments = voteSimEndogenousVoting_vary(parameters,rand_set)

    m = emp_moments .- simMoments
    y = m' * W * m

    return y
end

# ╔═╡ 1bb22460-37c4-11ec-3897-732e028fc413
begin

# Define the quasi-random starting guesses for the parameters. These are found in part B of the appendix

h0_v_in = rand(Uniform(0.2, 0.4), 1)
h0_nv_in = rand(Uniform(0.2, 0.4), 1)
r_v_in = rand(Uniform(0.2, 0.4), 1)
r_nv_in = rand(Uniform(0.2, 0.4), 1)
eta_v_in = rand(Uniform(0.0, 0.5), 1)
eta_nv_in = rand(Uniform(0.0, 0.5), 1)
mu_s_v_in = rand(Uniform(-50.0, 0.0), 1)
mu_s_nv_in = rand(Uniform(-50.0, 0.0), 1)
sigma_s_v_in = rand(Uniform(0.0, 50.0), 1)
sigma_s_nv_in = rand(Uniform(0.0, 50.0), 1)
S_svy_v_in = rand(Uniform(0.0, 10.0), 1)
S_svy_nv_in = rand(Uniform(0.0, 10.0), 1)
timeval_v_in = rand(Uniform(0.0, 100.0), 1)
timeval_nv_in = rand(Uniform(0.0, 100.0), 1)
mu_sv_in = rand(Uniform(-20.0, 20.0), 1)
mu_sn_in = rand(Uniform(-30.0, 10.0), 1)
sigma_svn_in = rand(Uniform(0.0, 30.0), 1)
L_in = rand(Uniform(0.0, 20.0), 1)
mu_eps_in = rand(Uniform(-30.0, 100.0), 1)
sigma_eps_in = rand(Uniform(50.0, 200.0), 1)

params_init=[h0_v_in,h0_nv_in,r_v_in,r_nv_in,eta_v_in,eta_nv_in,mu_s_v_in,mu_s_nv_in,sigma_s_v_in,
sigma_s_nv_in,S_svy_v_in,S_svy_nv_in,timeval_v_in,timeval_nv_in,mu_sv_in,mu_sn_in,sigma_svn_in,L_in,mu_eps_in,sigma_eps_in,]
params_init = [item for sublist in params_init for item in sublist] # flatten array

end;

# ╔═╡ e33132a0-37c5-11ec-1079-e3c9bba57502
begin

# !!! Read Before Running !!!

# This cell computes the estimates. On our machines it takes around 10-13 minutes for each 500 iterations. The number of iterations needed for convergence depends a lot on the starting guesses, nonetheless 5000 iterations are usually sufficient to get a wsse (weighted sum of squared errors) of 160 which is the one found also by the authors. In their Matlab code they run 720 optimization routines starting from different initial parameters to avoid being stuck in local minima and take as the best estimates those that provide the min wsse. Given time constraints we run the algorithm a dozen times and took the estimates that provided the lower weighted sum of squared errors

# !!! We reduced the number of iterations since Pluto would start this cell automatically and does not support stopping cells on windows. To find the estimates change the number of iterations from 10 to 5000 !!!

sol = optimize(x->criterion(x,rand_vec),params_init,NelderMead(),Optim.Options(
		       iterations = 10))
res = Optim.minimizer(sol)   # vector containing the solutions

end;

# ╔═╡ 9b7b0c70-38c8-11ec-1484-1b1f791d99dd
begin

# Since the estimation procedure can take a lot of time and depends on the rand_vec and the starting guesses, here we provide the data to replicate our best estimates. The file estimates2 contains the estimates found by the authors, our best estimates and the initial parameters we used to found them. random_vector_used contains the rand_vec we used to compute our best estimates

data, header = readdlm(raw"../input/estimates2.csv",',',header=true)
best_estimates = DataFrame(data, vec(header))

be = best_estimates[:,2]                # our best estimates

data, header = readdlm(raw"../input/random_vector_used.csv",',',header=true)
best_random_vec = DataFrame(data, vec(header))
rand_vec_used = [best_random_vec[:,1],best_random_vec[:,2],best_random_vec[:,3],best_random_vec[:,4]]
params_init_used = best_estimates[:,3]  # The random draws used to compute the best                                             estimates

end;

# ╔═╡ 5b42d010-38c9-11ec-1474-e103040023f8
begin

# This cell computes our best estimates given the rand_vec and the initial parameters used

# !!! change the number of iterations to 10.000 to find our best estimates !!!

sol = optimize(x->criterion(x, rand_vec_used), params_init_used, NelderMead(),
Optim.Options(iterations = 10)) # We set 10.000 iterations even if the algorithm will                                   converge before 10.000 iterations
be = Optim.minimizer(sol)

end

# ╔═╡ c8ce3502-38c6-11ec-3b68-f1513e5f03b5
begin

# Compare our best estimates with the ones obtained by the authors

# These are the authors' estimates
parameters_authors = [0.38,0.36,0.38,0.30,0.14,0.16,-22.6,-27.7,26.9,24.7,1.6,1.2,42.7,23.9,-3.9,-11.3,9.5,7.6,64.1,318.7]

# This difference is positive if the wsse obtained using our estimates is lower

dif = criterion(parameters_authors, rand_vec) - criterion(be, rand_vec)

end

# ╔═╡ 3f8c7950-38cb-11ec-1cf1-e79ce14eb776
# Below we compute the Jacobian using ForwardDiff but we do not use this jacobian to compute the standard errors since we get two columns of zeroes for the gradient with respect to mu_eps and sigma_eps and so the matrix is not invertible. This problem is either due to (1) mistakes in the function we wrote to compute the simulated moments or (2) the way in which julia computes the jacobian. Regarding the former possibility: using the same parameters and the same random vectors as the authors, we do see very small differences between ours and the authors' simulated moments (the largest difference being only 0.002); after inspecting it carefully, we cannot attribute these small discrepancies to mistakes in the function we wrote but it is possible that they are due to different default options between julia and matlab's operators (e.g., in rounding). Regarding the latter possibility: in julia, we computed the jacobian using the ForwardDiff Pkg while the authors used the Jacobianest function in Matlab; from a comment in the documentation of this Matlab function,  we learnt that "the error term on finite differences jacobian's estimates has a second order component, but also some 4th and 6th order terms in it" and the Jacobianest function uses Romberg extrapolation to improve the estimate, but ForwardDiff uses automatic differentiation so it should return a very precise jacobian.

jac_julia = ForwardDiff.jacobian(x->voteSimEndogenousVoting_vary(x, rand_vec_used), be);

# ╔═╡ 1f237b10-3e1f-11ec-2dde-092bc72714fb
# We load the jacobian evaluated in our best estimates computed in Matlab using the "Jacobianest" function (the same used by the authors).

jac_matlab = readdlm(raw"../input/jac_mat.csv",',', header=false);

# ╔═╡ 7a5da8fe-3b00-11ec-1f2c-7fc389373302
begin

# Compute the standard errors.

#For more informations on the formula for the standard errors we refer the reader to part B of the appendix or to the python notebook.

# We use the same notation as in their matlab code. sim_adjust is a scalar, DFDY_CSD_jacest is our jacobian (101x20 matrix) evaluated in the minimum, W is the weighting matrix (101x101 diagonal matrix), VCcontrol is the variance-covariance matrix of the empirical moments (a 101x101 matrix).

Jm_Js = 13197 / sim_voters # nr. empirical obs / nr. simulated obs
sim_adjust = 1 + Jm_Js
DFDY_CSD_jacest = jac_matlab
A = DFDY_CSD_jacest'*W*DFDY_CSD_jacest
B = DFDY_CSD_jacest'*W*(sim_adjust.*VCcontrol[1:101,1:101])*W*DFDY_CSD_jacest
VC = A\B/A
standard_errors = diag(VC).^0.5 # se are the square root of the main diagonal

end;

# ╔═╡ 427cce6e-37fa-11ec-2a6f-edc5c54a32aa
begin

# Create and save a csv file containing our estimates and their standard errors

# parameters
params = ["h0_v","h0_nv","r_v","r_nv","eta_v","eta_nv","mu_s_v","mu_s_nv","sigma_s_v","sigma_s_nv","S_svy_v","S_svy_nv","timeval_v","timeval_nv","mu_sv","mu_sn","sigma_si","L_in","mu_eps","sigma_eps"]

# authors' standard errors
se_authors= [0.0089,0.0092,0.0204,0.0183,0.1232,0.2987,2.9580,5.5938,5.5176,6.6687,1.2084,1.6149,9.4438,15.1539,1.4858,1.6422, 2.9750,1.5247,61.9195,248.6525]

table3=DataFrame()
table3.parameters = params
table3.our_point_est=round.(be, digits=2)
table3.our_se=round.(standard_errors, digits=4)
table3.authors_point_est=parameters_authors
table3.authors_se = se_authors

CSV.write(raw"../output/estimates_julia.csv", table3)
table3

end

# ╔═╡ Cell order:
# ╠═9c0d53a0-38c3-11ec-0d4d-1d549264f0f0
# ╠═92cfa5a2-5ce9-11ec-0e3e-61fe154ce557
# ╠═b34d84a0-5ce9-11ec-1316-39f597588bef
# ╠═73767a00-5db9-11ec-1bef-b9e6f6e711cc
# ╠═15e16bb0-3560-11ec-238e-e3018fb17cf9
# ╠═c6a68930-3588-11ec-08e1-7d0922b9816b
# ╠═c6731d20-3588-11ec-3127-51520038c40c
# ╠═c6468ee0-3588-11ec-0d41-cbfce972d21c
# ╠═c616cc50-3588-11ec-25bc-1b52b9c8694d
# ╠═1bb22460-37c4-11ec-3897-732e028fc413
# ╠═e33132a0-37c5-11ec-1079-e3c9bba57502
# ╠═9b7b0c70-38c8-11ec-1484-1b1f791d99dd
# ╠═5b42d010-38c9-11ec-1474-e103040023f8
# ╠═c8ce3502-38c6-11ec-3b68-f1513e5f03b5
# ╠═3f8c7950-38cb-11ec-1cf1-e79ce14eb776
# ╠═1f237b10-3e1f-11ec-2dde-092bc72714fb
# ╠═7a5da8fe-3b00-11ec-1f2c-7fc389373302
# ╠═427cce6e-37fa-11ec-2a6f-edc5c54a32aa

ext: jl