10  推定 1

using CSV, Tables
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

10.1 パラメタの設定

NumSimMarkets = 500;
NumSimPeriods = 50;

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
]...);

TransitionMat = [0.7 0.3; 0.4 0.6];

10.2 前章で生成したデータの読み込み

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

10.3 CCPの推定

EstimatedCCP1Mat, EstimatedCCP2Mat = estimateCCPMat(
    FakeData, global_param
);

10.4 遷移行列の推定

EstimatedTransition = estimateTransition(FakeData);

hcat(EstimatedCCP1Mat[:, 2], EstimatedCCP2Mat[:, 2])

8×2 Matrix{Float64}:
 0.671276  0.718293
 0.70041   0.800378
 0.827299  0.746308
 0.797659  0.767725
 0.705617  0.714142
 0.725334  0.761288
 0.789778  0.76567
 0.787342  0.736392

EstimatedTransition

2×2 Matrix{Float64}:
 0.700071  0.299929
 0.400289  0.599711

10.4.1 関数が正しく挙動しているかのチェック

output = updateCCP(
    parameterMat,
    CCP1Mat,
    CCP2Mat,
    TransitionMat,
    global_param
);


println("Difference between CCP in GDP and predicted CCP")

Difference between CCP in GDP and predicted CCP

CCP1Mat - output[1]

8×3 Matrix{Float64}:
  0.0         -1.95367e-9  1.95367e-9
  0.0         -2.42611e-9  2.42611e-9
 -1.33685e-9   1.33685e-9  0.0
 -1.75907e-9   1.75907e-9  0.0
  0.0         -1.94553e-9  1.94553e-9
  0.0         -2.46917e-9  2.46917e-9
 -1.3948e-9    1.3948e-9   0.0
 -2.09866e-9   2.09866e-9  0.0

CCP2Mat - output[2]

8×3 Matrix{Float64}:
  0.0         -6.91787e-9  6.91787e-9
 -5.08387e-9   5.08387e-9  0.0
  0.0         -5.43541e-9  5.43541e-9
 -5.30343e-9   5.30343e-9  0.0
  0.0         -6.37121e-9  6.37121e-9
 -5.46061e-9   5.46061e-9  0.0
  0.0         -5.44987e-9  5.44987e-9
 -6.33131e-9   6.33131e-9  0.0

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],
];

output_normalized = updateCCP(
    reshape(Normalized_TrueParam, (:, 2)),
    CCP1Mat,
    CCP2Mat,
    TransitionMat,
    global_param
);

println("Difference: predicted CCP in true parameter and normalized parameter")

Difference: predicted CCP in true parameter and normalized parameter

output_normalized[1] - output[1]

8×3 Matrix{Float64}:
 0.0           2.22045e-16  -1.66533e-16
 0.0           0.0           0.0
 2.498e-16    -2.22045e-16   0.0
 1.38778e-16  -1.11022e-16   0.0
 0.0           4.44089e-16  -3.33067e-16
 0.0           2.22045e-16  -2.77556e-16
 8.32667e-17  -2.22045e-16   0.0
 3.33067e-16  -2.22045e-16   0.0

output_normalized[2] - output[2]

8×3 Matrix{Float64}:
  0.0           1.11022e-16  -1.11022e-16
 -1.38778e-16   2.22045e-16   0.0
  0.0           0.0           0.0
  2.77556e-17   2.22045e-16   0.0
  0.0          -1.11022e-16   1.11022e-16
  0.0           0.0           0.0
  0.0           0.0          -1.66533e-16
 -5.55112e-17   1.11022e-16   0.0

10.5 Pesendorfer and Schmidt-Denglerの方法によるパラメタの推定

initial = [Parameters[1:4]; Parameters[6]];
mat_initial2 = initial .* reshape(rand(Uniform(0.6, 1.2), 50), (5, :));
result = zeros(6, 10);

@time for i in 1:10
    sol = optimize(
        x -> calculatePSDObjective(
            x,
            EstimatedCCP1Mat, 
            EstimatedCCP2Mat, 
            EstimatedTransition, 
            global_param
            ),
        mat_initial2[:, i],
        Optim.Options(show_trace = false)
    );
    result[:, i] = [sol.minimizer; sol.minimum];
end

result_pick = result[1:5, argmin(result[6, :])];

  0.894520 seconds (5.11 M allocations: 392.116 MiB, 8.88% gc time, 81.60% compilation time)

normalized_trueparam = [
    Normalized_TrueParam[1],
    Normalized_TrueParam[6],
    Normalized_TrueParam[2],
    Normalized_TrueParam[3],
    Normalized_TrueParam[5],
];

DataFrame(
    param_est = result_pick,
    normalized_trueparam = normalized_trueparam
)

5×2 DataFrame

Row	param_est	normalized_trueparam
	Float64	Float64
1	0.348787	0.3375
2	0.223868	0.2375
3	-0.291698	-0.27
4	0.492655	0.45
5	-2.22346	-2.25

numBootSample = 100;

bootindex = reshape(
    sample(1:NumSimMarkets, NumSimMarkets * numBootSample),
    (NumSimMarkets, 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 = zeros(NumSimPeriods * NumSimMarkets, 8);
    for m in 1:NumSimMarkets
        temp = FakeData[FakeData[:, 1] .== bootindex[m, b], :];
        bootsample[(1 + NumSimPeriods * (m - 1)):(NumSimPeriods * m), :] = temp;
    end

    output = Estimation_PS_bootstrap(bootsample, global_param);

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

  3.938225 seconds (17.45 M allocations: 13.188 GiB, 16.14% gc time, 1.48% compilation time)

println("CCP for firm 1")
DataFrame(
    "True values" => CCP1Mat[:, 2],
    "Estimated" => vec(EstimatedCCP1Mat[:, 2]),
    "SE" => vec(std(bootresult_CCP1, dims=2))
)

CCP for firm 1

8×3 DataFrame

Row	True values	Estimated	SE
	Float64	Float64	Float64
1	0.683563	0.671276	0.00916724
2	0.712662	0.70041	0.00767212
3	0.815651	0.827299	0.00620892
4	0.793819	0.797659	0.00577186
5	0.708067	0.705617	0.00959225
6	0.735965	0.725334	0.00899726
7	0.797547	0.789778	0.00860466
8	0.774137	0.787342	0.00711149

println("CCP for firm 2")
DataFrame(
    "True values" => CCP2Mat[:, 2],
    "Estimated" => vec(EstimatedCCP2Mat[:, 2]),
    "SE" => vec(std(bootresult_CCP2, dims=2))
)

CCP for firm 2

8×3 DataFrame

Row	True values	Estimated	SE
	Float64	Float64	Float64
1	0.714725	0.718293	0.00968694
2	0.792409	0.800378	0.00731078
3	0.741697	0.746308	0.00572638
4	0.768731	0.767725	0.00640628
5	0.737685	0.714142	0.0104247
6	0.772706	0.761288	0.00832577
7	0.76353	0.76567	0.00721454
8	0.7473	0.736392	0.00769535

println("Transition Probability (GG and BB)")
DataFrame(
    "True values" => diag(TransitionMat),
    "Estimated" => diag(EstimatedTransition),
    "SE" => vec(std(bootresult_transition, dims=2))
)

Transition Probability (GG and BB)

2×3 DataFrame

Row	True values	Estimated	SE
	Float64	Float64	Float64
1	0.7	0.700071	0.00337154
2	0.6	0.599711	0.0045105

DataFrame(
    "True values" => initial,
    "Normalized true values" => Normalized_TrueParam[[1, 6, 2, 3, 5]],
    "Estimated" => result_pick,
    "SE" => vec(std(bootresult_payoff, dims=2))
    )

5×4 DataFrame

Row	True values	Normalized true values	Estimated	SE
	Float64	Float64	Float64	Float64
1	0.3	0.3375	0.348787	0.0544961
2	0.2	0.2375	0.223868	0.0557081
3	-0.27	-0.27	-0.291698	0.0385321
4	0.45	0.45	0.492655	0.089435
5	-2.1	-2.25	-2.22346	0.0196745

@save "tmp/dynamic_game/estimates_1.jld2" EstimatedCCP1Mat EstimatedCCP2Mat EstimatedTransition;