11  推定 2

using DataFrames
using Random
using Distributions
using JLD2
using Statistics
using StatsBase
using LinearAlgebra
using Optim

mutable struct global_param_struct
    beta::Float64
    eulergamma::Float64
    entryMat::Matrix{Int}
    numFirmsVec::Vector{Int}
    economyVec::Vector{Int}
    possibleActionArray::Array{Int, 3}
end

for file in readdir("functions/dynamic_game/")
    include("functions/dynamic_game/" * file)
end

11.1 パラメタの設定

TransitionMat = [0.7 0.3; 0.4 0.6]

Random.seed!(123);

Parameters = [
    0.3,
    0.2,
    -0.27,
    0.45,
    -0.15,
    -2.10
];

parameterMat = hcat([
    [Parameters[i]; Parameters[3:6]]
    for i in 1:2
]...);

NumSimPeriods = 100;
NumSimFirms = 2;
NumSimulations = 1000;
NumSimMarkets = 8;

NumSimMarketsInData = 500;

numBootSample = 100;

11.2 前章で生成したデータと推定したパラメタの読み込み

@load "tmp/dynamic_game/data_workspace.jld2";
@load "tmp/dynamic_game/estimates_1.jld2";

InitialState = hcat(1:8, 1:8);;

EVrandom = reshape(
    -log.(-log.(rand(NumSimMarkets * NumSimPeriods * NumSimFirms * NumSimulations * 8 * 3))),
    (NumSimMarkets, NumSimPeriods, NumSimFirms, NumSimulations, 8, 3)
);

UNIrandom = reshape(
    rand(NumSimMarkets * NumSimPeriods * NumSimulations),
    (NumSimMarkets, NumSimPeriods, NumSimulations)
);

@time Wstar = VSigmaGeneration(
    CCP1Mat[:, 2],
    CCP2Mat[:, 2],
    TransitionMat,
    EVrandom,
    UNIrandom,
    InitialState,
    NumSimMarkets,
    NumSimulations,
    NumSimPeriods,
    global_param
);

W1star = Wstar[1];
W2star = Wstar[2];

  1.297144 seconds (6.04 M allocations: 569.625 MiB, 9.94% gc time, 59.58% compilation time)

param1 = [parameterMat[:, 1]; 1];
param2 = [parameterMat[:, 2]; 1];

DataFrame(
    "True value" => exanteV1,
    "Simulated value" => W1star' * param1,
    "Difference" => exanteV1 - W1star' * param1
)

8×3 DataFrame

Row	True value	Simulated value	Difference
	Float64	Float64	Float64
1	4.6568	4.53718	0.119622
2	4.58732	4.63857	-0.0512471
3	6.56876	6.52198	0.0467797
4	6.15887	6.18534	-0.0264716
5	4.61238	4.60842	0.00396147
6	4.54617	4.58591	-0.0397407
7	6.01581	5.87433	0.141481
8	5.6088	5.69862	-0.0898211

DataFrame(
    "True value" => exanteV2,
    "Simulated value" => W2star' * param2,
    "Difference" => exanteV2 - W2star' * param2
)

8×3 DataFrame

Row	True value	Simulated value	Difference
	Float64	Float64	Float64
1	4.44119	4.48669	-0.0454969
2	6.17892	6.14646	0.0324623
3	4.38004	4.47488	-0.0948405
4	5.77645	5.80853	-0.03208
5	4.40139	4.47127	-0.0698898
6	5.62999	5.54035	0.0896341
7	4.34285	4.40536	-0.0625095
8	5.22972	5.20872	0.0209918

Normalized_TrueParam = [
    Parameters[1] - ((1 - global_param.beta) / global_param.beta) * Parameters[5],
    Parameters[3],
    Parameters[4],
    0,
    Parameters[6] + Parameters[5],
    Parameters[2] - ((1 - global_param.beta) / global_param.beta) * Parameters[5],
    Parameters[3],
    Parameters[4],
    0,
    Parameters[6] + Parameters[5],
];

normparam1 = [Normalized_TrueParam[1:5]; 1];
normparam2 = [Normalized_TrueParam[6:10]; 1];

DataFrame(
    "True value" => exanteV1,
    "Simulated value (original param)" => W1star' * param1,
    "Simulated value (normalized param)" => W1star' * normparam1
)

8×3 DataFrame

Row	True value	Simulated value (original param)	Simulated value (normalized param)
	Float64	Float64	Float64
1	4.6568	4.53718	4.53718
2	4.58732	4.63857	4.63857
3	6.56876	6.52198	6.70948
4	6.15887	6.18534	6.37284
5	4.61238	4.60842	4.60842
6	4.54617	4.58591	4.58591
7	6.01581	5.87433	6.06183
8	5.6088	5.69862	5.88612

DataFrame(
    "True value" => exanteV2,
    "Simulated value (original param)" => W2star' * param2,
    "Simulated value (normalized param)" => W2star' * normparam2
)

8×3 DataFrame

Row	True value	Simulated value (original param)	Simulated value (normalized param)
	Float64	Float64	Float64
1	4.44119	4.48669	4.48669
2	6.17892	6.14646	6.33396
3	4.38004	4.47488	4.47488
4	5.77645	5.80853	5.99603
5	4.40139	4.47127	4.47127
6	5.62999	5.54035	5.72785
7	4.34285	4.40536	4.40536
8	5.22972	5.20872	5.39622

println("Difference b/w two simulated values (with and without normalization)")
DataFrame(
    "Firm 1" => W1star' * param1 - W1star' * normparam1,
    "Firm 2" => W2star' * param2 - W2star' * normparam2,
)

Difference b/w two simulated values (with and without normalization)

8×2 DataFrame

Row	Firm 1	Firm 2
	Float64	Float64
1	2.15419e-11	1.96705e-11
2	2.2764e-11	-0.1875
3	-0.1875	2.10454e-11
4	-0.1875	-0.1875
5	2.31068e-11	1.90594e-11
6	2.28786e-11	-0.1875
7	-0.1875	2.11973e-11
8	-0.1875	-0.1875

11.3 Forward simulationを用いたPesendorfer and Schmidt-Denglerによるパラメタの推定

@time Wstar = VSigmaGeneration(
    EstimatedCCP1Mat[:, 2],
    EstimatedCCP2Mat[:, 2],
    EstimatedTransition,
    EVrandom,
    UNIrandom,
    InitialState,
    NumSimMarkets,
    NumSimulations,
    NumSimPeriods,
    global_param
);

W1star = Wstar[1];
W2star = Wstar[2];

  0.278932 seconds (4.80 M allocations: 489.018 MiB, 19.05% gc time)

initial = [0.3375, 0.2375, -0.27, 0.45, -2.25];

@time resultForwardPSD = optimize(
    x -> Estimation_forward_PSD(
        x,
        W1star,
        W2star,
        EstimatedTransition,
        EstimatedCCP1Mat,
        EstimatedCCP2Mat,
        global_param
    ),
    initial,
    Optim.Options(show_trace = false)
)

bootindex = reshape(
    sample(1:NumSimMarketsInData, NumSimMarketsInData * numBootSample),
    (NumSimMarketsInData, numBootSample)
);

bootresult_transition = zeros(2, numBootSample);
bootresult_CCP1 = zeros(8, numBootSample);
bootresult_CCP2 = zeros(8, numBootSample);
bootresult_payoff = zeros(5, numBootSample);

@time for b in 1:numBootSample

    bootsample = vcat(
        [FakeData[FakeData[:, 1] .== bootindex[m, b], :] for m in 1:NumSimMarketsInData]...
        );

    output = Bootstrap_PS_forward(
        bootsample, 
        global_param,
        EVrandom,
        UNIrandom,
        InitialState,
        NumSimMarkets,
        NumSimulations,
        NumSimPeriods
        );

    bootresult_payoff[:, b] = output[1].minimizer;
    bootresult_CCP1[:, b] = output[2][:, 2]
    bootresult_CCP2[:, b] = output[3][:, 2]
    bootresult_transition[:, b] = diag(output[4]);

end

DataFrame(
    "True param" => [Parameters[1:4]; Parameters[6]],
    "True param (normalized)" => Normalized_TrueParam[[1, 6, 2, 3, 5]],
    "Estimates" => resultForwardPSD.minimizer,
    "SE" => vec(std(bootresult_payoff, dims = 2))
)

5×4 DataFrame

Row	True param	True param (normalized)	Estimates	SE
	Float64	Float64	Float64	Float64
1	0.3	0.3375	0.394852	0.0544667
2	0.2	0.2375	0.319143	0.0561027
3	-0.27	-0.27	-0.319576	0.0373488
4	0.45	0.45	0.440578	0.0926289
5	-2.1	-2.25	-2.22364	0.0186967

11.4 BBLの不等式推定量を用いたパラメタの推定

NumPerturbations = 200;

PerturbedCCP1 = reshape(
    clamp.(
        repeat(EstimatedCCP1Mat[:, 2], NumPerturbations) + 
        rand(Normal(0, .1), 8 * NumPerturbations),
        0.001, 0.999
    ),
    (8, NumPerturbations)
);

PerturbedCCP2 = reshape(
    clamp.(
        repeat(EstimatedCCP2Mat[:, 2], NumPerturbations) + 
        rand(Normal(0, .1), 8 * NumPerturbations),
        0.001, 0.999
    ),
    (8, NumPerturbations)
);

W1_all = zeros(6, NumSimMarkets, NumPerturbations);
W2_all = zeros(6, NumSimMarkets, NumPerturbations);

@time for per in 1:NumPerturbations

    W1_p = VSigmaGeneration(
        PerturbedCCP1[:, per],
        EstimatedCCP2Mat[:, 2],
        EstimatedTransition,
        EVrandom,
        UNIrandom,
        InitialState,
        NumSimMarkets,
        NumSimulations,
        NumSimPeriods,
        global_param
    );
    W1_all[:, :, per] = W1_p[1];

    W2_p = VSigmaGeneration(
        EstimatedCCP1Mat[:, 2],
        PerturbedCCP2[:, per],
        EstimatedTransition,
        EVrandom,
        UNIrandom,
        InitialState,
        NumSimMarkets,
        NumSimulations,
        NumSimPeriods,
        global_param
    );
    W2_all[:, :, per] = W2_p[2];

end

@time BBLobjective(
        initial,
        NumPerturbations,
        W1star,
        W2star,
        W1_all,
        W2_all
    )

initial = [0.3, 0.2, -0.27, 0.45, -2.1];

@time resultBBL = optimize(
    x -> BBLobjective(
        x,
        NumPerturbations,
        W1star,
        W2star,
        W1_all,
        W2_all
    ),
    initial,
    Optim.Options(show_trace = false)
)

bootindex = reshape(
    sample(1:NumSimMarketsInData, NumSimMarketsInData * numBootSample),
    (NumSimMarketsInData, numBootSample)
);

bootresult_transition_BBL = zeros(2, numBootSample);
bootresult_CCP1_BBL = zeros(8, numBootSample);
bootresult_CCP2_BBL = zeros(8, numBootSample);
bootresult_payoff_BBL = zeros(5, numBootSample);

@time for b in 1:numBootSample

    bootsample = vcat(
        [FakeData[FakeData[:, 1] .== bootindex[m, b], :] for m in 1:NumSimMarketsInData]...
        );

    output = Bootstrap_BBL(
        bootsample, 
        global_param,
        EVrandom,
        UNIrandom,
        InitialState,
        NumSimMarkets,
        NumSimulations,
        NumSimPeriods,
        NumPerturbations
        );

    bootresult_payoff_BBL[:, b] = output[1].minimizer;
    bootresult_CCP1_BBL[:, b] = output[2][:, 2]
    bootresult_CCP2_BBL[:, b] = output[3][:, 2]
    bootresult_transition_BBL[:, b] = diag(output[4]);

end

DataFrame(
    "True param" => [Parameters[1:4]; Parameters[6]],
    "True param (normalized)" => Normalized_TrueParam[[1, 6, 2, 3, 5]],
    "Estimates" => resultBBL.minimizer,
    "SE" => vec(std(bootresult_payoff_BBL, dims = 2))
)

5×4 DataFrame

Row	True param	True param (normalized)	Estimates	SE
	Float64	Float64	Float64	Float64
1	0.3	0.3375	0.361364	0.156315
2	0.2	0.2375	0.170143	0.167567
3	-0.27	-0.27	-0.292871	0.151854
4	0.45	0.45	0.525381	0.262609
5	-2.1	-2.25	-2.10009	0.0884199

DataFrame(
    "True param" => [Parameters[1:4]; Parameters[6]],
    "True param (normalized)" => Normalized_TrueParam[[1, 6, 2, 3, 5]],
    "Estimates (Forward P-SD)" => resultForwardPSD.minimizer,
    "SE (Forward P-SD)" => vec(std(bootresult_payoff, dims = 2)),
    "Estimates (BBL)" => resultBBL.minimizer,
    "SE (BBL)" => vec(std(bootresult_payoff_BBL, dims = 2))
)

5×6 DataFrame

Row	True param	True param (normalized)	Estimates (Forward P-SD)	SE (Forward P-SD)	Estimates (BBL)	SE (BBL)
	Float64	Float64	Float64	Float64	Float64	Float64
1	0.3	0.3375	0.394852	0.0544667	0.361364	0.156315
2	0.2	0.2375	0.319143	0.0561027	0.170143	0.167567
3	-0.27	-0.27	-0.319576	0.0373488	-0.292871	0.151854
4	0.45	0.45	0.440578	0.0926289	0.525381	0.262609
5	-2.1	-2.25	-2.22364	0.0186967	-2.10009	0.0884199