9  疑似データの生成

using Random
using StatsBase
using CSV, Tables
using LinearAlgebra
using JLD2

mutable struct array_struct
    CCP1Mat::Matrix{Float64}
    CCP2Mat::Matrix{Float64}
    stateTransitionMat::Matrix{Float64}
    equilibriumProfitFirm1::Matrix{Float64}
    equilibriumProfitFirm2::Matrix{Float64}
    expectedShockUnderBestActionFirm1::Matrix{Float64}
    expectedShockUnderBestActionFirm2::Matrix{Float64}
    exanteV1::Vector{Float64}
    exanteV2::Vector{Float64}
    updatedCCP1Numerator::Matrix{Float64}
    updatedCCP2Numerator::Matrix{Float64}
    CCP1MatUpdated::Matrix{Float64}
    CCP2MatUpdated::Matrix{Float64}
    stateTransitionMatFirm2ActionSpecific::Array{Float64, 3}
    stateTransitionMatFirm1ActionSpecific::Array{Float64, 3}
    exanteV1Updated::Vector{Float64}
    exanteV2Updated::Vector{Float64}
end

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

9.1 パラメタの設定

Random.seed!(123);

beta = 0.8;
eulergamma = Float64(MathConstants.eulergamma);
TransitionMat = [0.7 0.3; 0.4 0.6]

entryMat = repeat([
    0 0;
    0 1;
    1 0;
    1 1;
], 2);
numFirmsVec = vec(sum(entryMat, dims = 2));
economyVec = [repeat([1], 4); repeat([0], 4)];

possibleActionArray = zeros(Int, 8, 3, 2);
possibleActionArray[:, :, 1] .= repeat([
    0 1 1;
    0 1 1;
    1 1 0;
    1 1 0
], 2);
possibleActionArray[:, :, 2] .= repeat([
    0 1 1;
    1 1 0;
    0 1 1;
    1 1 0
], 2);

global_param = global_param_struct(
    beta,
    eulergamma,
    entryMat,
    numFirmsVec,
    economyVec,
    possibleActionArray
);

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

9.2 各企業の利潤の計算

profitFirm1 = calculateFirmProfit(
    parameterMat,
    global_param.entryMat,
    global_param.numFirmsVec,
    global_param.economyVec,
    global_param.possibleActionArray,
    1
);
profitFirm2 = calculateFirmProfit(
    parameterMat,
    global_param.entryMat,
    global_param.numFirmsVec,
    global_param.economyVec,
    global_param.possibleActionArray,
    2
);

9.3 CCPの計算

CCP1Stay = repeat([0.5], 8);
CCP2Stay = repeat([0.5], 8);

CCP1Mat = (
    hcat(
        1 .- CCP1Stay,
        CCP1Stay,
        1 .- CCP1Stay,
    ) .* 
    global_param.possibleActionArray[:, :, 1]
    );
CCP2Mat = (
    hcat(
        1 .- CCP2Stay,
        CCP2Stay,
        1 .- CCP2Stay,
    ) .* 
    global_param.possibleActionArray[:, :, 2]
    );

stateTransitionMat = generateStateTransitionMat(
    TransitionMat, 
    CCP1Mat, 
    CCP2Mat, 
    global_param.possibleActionArray
    );

equilibriumProfitFirm1 = calculateProfitGivenOpponentCCP(profitFirm1, CCP2Mat);
equilibriumProfitFirm2 = calculateProfitGivenOpponentCCP(profitFirm2, CCP1Mat);

expectedShockUnderBestActionFirm1 = (
    eulergamma .- replace(log.(CCP1Mat), -Inf => 0)
)
expectedShockUnderBestActionFirm2 = (
    eulergamma .- replace(log.(CCP2Mat), -Inf => 0)
)

exanteV1 = vec(
    (I(8) - beta .* stateTransitionMat) \ 
    sum(CCP1Mat .* (equilibriumProfitFirm1 + expectedShockUnderBestActionFirm1), dims = 2)
);

exanteV2 = vec(
    (I(8) - beta .* stateTransitionMat) \ 
    sum(CCP2Mat .* (equilibriumProfitFirm2 + expectedShockUnderBestActionFirm2), dims = 2)
);

# updatedCCP1Numerator = zeros(8, 3);
# updatedCCP2Numerator = zeros(8, 3);
# CCP1MatUpdated = zeros(8, 3);
# CCP2MatUpdated = zeros(8, 3);
# stateTransitionMatFirm2ActionSpecific = zeros(8, 8, 3);
# stateTransitionMatFirm1ActionSpecific = zeros(8, 8, 3);

data_generate_struct = array_struct(
    CCP1Mat,
    CCP2Mat,
    stateTransitionMat,
    equilibriumProfitFirm1,
    equilibriumProfitFirm2,
    expectedShockUnderBestActionFirm1,
    expectedShockUnderBestActionFirm2,
    exanteV1,
    exanteV2,
    zeros(8, 3),
    zeros(8, 3),
    zeros(8, 3),
    zeros(8, 3),
    zeros(8, 8, 3),
    zeros(8, 8, 3),
    zeros(8),
    zeros(8)
);

diffExanteV = 1;
iter = 0;

@time while (diffExanteV > 1e-12)

    calculateFirm2ActionSpecificTransitionMat!(
        TransitionMat, 
        data_generate_struct.CCP2Mat,
        data_generate_struct.stateTransitionMatFirm2ActionSpecific,
        global_param.possibleActionArray
        );
    calculateFirm1ActionSpecificTransitionMat!(
        TransitionMat, 
        data_generate_struct.CCP1Mat,
        data_generate_struct.stateTransitionMatFirm1ActionSpecific,
        global_param.possibleActionArray
        );

    for i in 1:3
        data_generate_struct.updatedCCP1Numerator[:, i] .= exp.(
            view(data_generate_struct.equilibriumProfitFirm1, :, i) .+ 
            beta .* view(data_generate_struct.stateTransitionMatFirm2ActionSpecific, :, :, i) * 
            data_generate_struct.exanteV1 .*
            view(global_param.possibleActionArray, :, i, 1)
        );
    end
    data_generate_struct.CCP1MatUpdated .= (
        (
            data_generate_struct.updatedCCP1Numerator ./ 
            (sum(data_generate_struct.updatedCCP1Numerator, dims = 2) .- 1)) .* 
        view(global_param.possibleActionArray, :, :, 1)
    );

    for i in 1:3
        data_generate_struct.updatedCCP2Numerator[:, i] .= exp.(
            view(data_generate_struct.equilibriumProfitFirm2, :, i) .+ 
            beta .* view(data_generate_struct.stateTransitionMatFirm1ActionSpecific, :, :, i) * 
            data_generate_struct.exanteV2 .*
            view(global_param.possibleActionArray, :, i, 2)
        );
    end
    data_generate_struct.CCP2MatUpdated .= (
        (
            data_generate_struct.updatedCCP2Numerator ./ 
            (sum(data_generate_struct.updatedCCP2Numerator, dims = 2) .- 1)
            ) .* 
        view(global_param.possibleActionArray, :, :, 2)
    );

    data_generate_struct.stateTransitionMat .= generateStateTransitionMat(
        TransitionMat, 
        data_generate_struct.CCP1MatUpdated, 
        data_generate_struct.CCP2MatUpdated, 
        global_param.possibleActionArray
        );

    data_generate_struct.equilibriumProfitFirm1 .= calculateProfitGivenOpponentCCP(
        profitFirm1, 
        data_generate_struct.CCP2MatUpdated
        );
    data_generate_struct.equilibriumProfitFirm2 .= calculateProfitGivenOpponentCCP(
        profitFirm2, 
        data_generate_struct.CCP1MatUpdated
        );

    data_generate_struct.expectedShockUnderBestActionFirm1 .= (
        eulergamma .- replace(log.(data_generate_struct.CCP1MatUpdated), -Inf => 0)
    )
    data_generate_struct.expectedShockUnderBestActionFirm2 .= (
        eulergamma .- replace(log.(data_generate_struct.CCP2MatUpdated), -Inf => 0)
    )

    data_generate_struct.exanteV1Updated .= vec(
        (I(8) - beta .* data_generate_struct.stateTransitionMat) \ 
        sum(
            data_generate_struct.CCP1MatUpdated .* (
                data_generate_struct.equilibriumProfitFirm1 + 
                data_generate_struct.expectedShockUnderBestActionFirm1
                ), dims = 2)
    );

    data_generate_struct.exanteV2Updated .= vec(
        (I(8) - beta .* data_generate_struct.stateTransitionMat) \ 
        sum(
            data_generate_struct.CCP2MatUpdated .* (
                data_generate_struct.equilibriumProfitFirm2 + 
                data_generate_struct.expectedShockUnderBestActionFirm2
                ), dims = 2)
    );

    diffExanteV = sum(
        (data_generate_struct.exanteV1Updated - data_generate_struct.exanteV1).^2
        ) + sum(
            (data_generate_struct.exanteV2Updated - data_generate_struct.exanteV2).^2
            );

    data_generate_struct.exanteV1 .= data_generate_struct.exanteV1Updated[:];
    data_generate_struct.exanteV2 .= data_generate_struct.exanteV2Updated[:];

    data_generate_struct.CCP1Mat .= data_generate_struct.CCP1MatUpdated[:, :];
    data_generate_struct.CCP2Mat .= data_generate_struct.CCP2MatUpdated[:, :];

    println(diffExanteV)
    iter += 1;

end

CCP1Mat .= data_generate_struct.CCP1Mat;
CCP2Mat .= data_generate_struct.CCP2Mat;

exanteV1 = data_generate_struct.exanteV1;
exanteV2 = data_generate_struct.exanteV2;

9.059385123406884
0.0010509979269102616
5.0928397934340845e-6
5.6083848572020334e-8
5.889460824473882e-10
5.5728209722870154e-12
3.866293900916097e-14
  1.262130 seconds (5.23 M allocations: 332.698 MiB, 3.42% gc time, 98.65% compilation time)

9.4 疑似データの生成

NumSimMarkets = 500;
NumSimPeriods = 50;
NumSimFirms = 2;

InitialState = sample(1:8, NumSimMarkets);

RandomNumbers = reshape(
    rand(NumSimMarkets * NumSimPeriods * (NumSimFirms + 1)),
    (NumSimMarkets, NumSimPeriods, NumSimFirms + 1)
);

FakeData = zeros(NumSimMarkets * NumSimPeriods, 8);

for m in 1:NumSimMarkets, t in 1:NumSimPeriods

    FakeData[(m - 1) * NumSimPeriods + t, 1] = m;
    FakeData[(m - 1) * NumSimPeriods + t, 2] = t;

    if t == 1

        FakeData[(m - 1) * NumSimPeriods + 1, 3] = InitialState[m];

        FakeData[(m - 1) * NumSimPeriods + 1, 4] = (
            (InitialState[m] >= 1) & (InitialState[m] <= 4) ? 1 :
            2
        );

    else

        sprev = FakeData[(m - 1) * NumSimPeriods + t - 1, 3]
        a1prev = FakeData[(m - 1) * NumSimPeriods + t - 1, 7]
        a2prev = FakeData[(m - 1) * NumSimPeriods + t - 1, 8]

        if (sprev >= 1) & (sprev <= 4)
            FakeData[(m - 1) * NumSimPeriods + t, 4] = (
                (RandomNumbers[m, t, 3] < TransitionMat[1, 1]) ? 1 :
                2
            );
        else
            FakeData[(m - 1) * NumSimPeriods + t, 4] = (
                (RandomNumbers[m, t, 3] < TransitionMat[2, 2]) ? 2 :
                1
            );
        end

        snow = FakeData[(m - 1) * NumSimPeriods + t, 4];
        n1t = FakeData[(m - 1) * NumSimPeriods + t - 1, 5] + a1prev;
        n2t = FakeData[(m - 1) * NumSimPeriods + t - 1, 6] + a2prev;

        FakeData[(m - 1) * NumSimPeriods + t, 3] = (
            (snow == 1) & (n1t == 0) & (n2t == 0) ? 1 :
            (snow == 1) & (n1t == 0) & (n2t == 1) ? 2 :
            (snow == 1) & (n1t == 1) & (n2t == 0) ? 3 :
            (snow == 1) & (n1t == 1) & (n2t == 1) ? 4 :
            (snow == 2) & (n1t == 0) & (n2t == 0) ? 5 :
            (snow == 2) & (n1t == 0) & (n2t == 1) ? 6 :
            (snow == 2) & (n1t == 1) & (n2t == 0) ? 7 :
            (snow == 2) & (n1t == 1) & (n2t == 1) ? 8 :
            0
        );

    end

    s = Int(FakeData[(m - 1) * NumSimPeriods + t, 3]);

    if s == 1
        FakeData[(m - 1) * NumSimPeriods + t, 4:6] .= [1, 0, 0];
        if (RandomNumbers[m, t, 1] > CCP1Mat[s, 2])
            FakeData[(m - 1) * NumSimPeriods + t, 7] = 1;
        end
        if (RandomNumbers[m, t, 2] > CCP2Mat[s, 2])
            FakeData[(m - 1) * NumSimPeriods + t, 8] = 1;
        end
    elseif s == 2
        FakeData[(m - 1) * NumSimPeriods + t, 4:6] .= [1, 0, 1];
        if (RandomNumbers[m, t, 1] > CCP1Mat[s, 2])
            FakeData[(m - 1) * NumSimPeriods + t, 7] = 1;
        end
        if (RandomNumbers[m, t, 2] > CCP2Mat[s, 2])
            FakeData[(m - 1) * NumSimPeriods + t, 8] = -1;
        end
    elseif s == 3
        FakeData[(m - 1) * NumSimPeriods + t, 4:6] .= [1, 1, 0];
        if (RandomNumbers[m, t, 1] > CCP1Mat[s, 2])
            FakeData[(m - 1) * NumSimPeriods + t, 7] = -1;
        end
        if (RandomNumbers[m, t, 2] > CCP2Mat[s, 2])
            FakeData[(m - 1) * NumSimPeriods + t, 8] = 1;
        end
    elseif s == 4
        FakeData[(m - 1) * NumSimPeriods + t, 4:6] .= [1, 1, 1];
        if (RandomNumbers[m, t, 1] > CCP1Mat[s, 2])
            FakeData[(m - 1) * NumSimPeriods + t, 7] = -1;
        end
        if (RandomNumbers[m, t, 2] > CCP2Mat[s, 2])
            FakeData[(m - 1) * NumSimPeriods + t, 8] = -1;
        end
    elseif s == 5
        FakeData[(m - 1) * NumSimPeriods + t, 4:6] .= [2, 0, 0];
        if (RandomNumbers[m, t, 1] > CCP1Mat[s, 2])
            FakeData[(m - 1) * NumSimPeriods + t, 7] = 1;
        end
        if (RandomNumbers[m, t, 2] > CCP2Mat[s, 2])
            FakeData[(m - 1) * NumSimPeriods + t, 8] = 1;
        end
    elseif s == 6
        FakeData[(m - 1) * NumSimPeriods + t, 4:6] .= [2, 0, 1];
        if (RandomNumbers[m, t, 1] > CCP1Mat[s, 2])
            FakeData[(m - 1) * NumSimPeriods + t, 7] = 1;
        end
        if (RandomNumbers[m, t, 2] > CCP2Mat[s, 2])
            FakeData[(m - 1) * NumSimPeriods + t, 8] = -1;
        end
    elseif s == 7
        FakeData[(m - 1) * NumSimPeriods + t, 4:6] .= [2, 1, 0];
        if (RandomNumbers[m, t, 1] > CCP1Mat[s, 2])
            FakeData[(m - 1) * NumSimPeriods + t, 7] = -1;
        end
        if (RandomNumbers[m, t, 2] > CCP2Mat[s, 2])
            FakeData[(m - 1) * NumSimPeriods + t, 8] = 1;
        end
    elseif s == 8
        FakeData[(m - 1) * NumSimPeriods + t, 4:6] .= [2, 1, 1];
        if (RandomNumbers[m, t, 1] > CCP1Mat[s, 2])
            FakeData[(m - 1) * NumSimPeriods + t, 7] = -1;
        end
        if (RandomNumbers[m, t, 2] > CCP2Mat[s, 2])
            FakeData[(m - 1) * NumSimPeriods + t, 8] = -1;
        end
    end
end

jldsave(
    "tmp/dynamic_game/data_workspace.jld2" ;
    FakeData, CCP1Mat, CCP2Mat, 
    exanteV1, exanteV2, global_param
)