5  基礎編

using CSV
using DataFrames
using StringEncodings
using LinearAlgebra
using Optim
using ForwardDiff
using Distributions

for file in ["calculateBRObjective.jl"]
    include("functions/entry_exit/" * file)
end

5.1 データの読み込み

data = CSV.File(
    open(read, "data/entry_exit_basic/MRIData.csv", enc"UTF-8"),
    missingstring = ["NA", ""],
    ) |> DataFrame
first(data, 5)

5×16 DataFrame

Row	HospitalID	CityCode	Kyukyu	Kinou	Sien	Hyoka	DepNeurology	DepNeurosurgery	NumBeds	MRINumOwn	Tesla	Population	Menseki	PopDensity	TaxableIncome	MRIOwnDum
	Int64	Int64	Int64?	Int64?	Int64?	Int64?	Int64?	Int64?	Int64?	Int64	Int64?	Int64?	Float64?	Float64?	Int64?	Int64
1	1	1101	1	0	0	0	0	0	243	1	2	220189	46.42	4743.4	3017	1
2	2	1101	0	0	0	0	0	0	180	1	0	220189	46.42	4743.4	3017	1
3	3	1101	1	0	0	0	0	0	110	1	1	220189	46.42	4743.4	3017	1
4	4	1101	1	0	0	0	1	1	134	3	2	220189	46.42	4743.4	3017	1
5	5	1101	0	0	0	0	0	0	250	1	2	220189	46.42	4743.4	3017	1

dataset = combine(
    groupby(data, :CityCode),
    :MRIOwnDum => sum => :NumMRI,
    [:Population, :Menseki, :PopDensity, :TaxableIncome] .=> mean .=> [:Pop, :Menseki, :PopDen, :Income],   
);
dataset[!, :Pop] = dataset.Pop / 1e+6;
dropmissing!(dataset, :);
numObs = nrow(dataset)

1459

5.2 Bresnahan and Reiss (1991) モデルの推定

hospitalNumCap = 6;
dataset[dataset.NumMRI .> hospitalNumCap, :NumMRI] .= hospitalNumCap;

obj_for_Optim = TwiceDifferentiable(
    x -> calculateBRObjective(x, dataset, hospitalNumCap),
    ones(hospitalNumCap + 1);
    autodiff = :forward
);
@time optim_res = optimize(
    obj_for_Optim,
    zeros(hospitalNumCap + 1),
    repeat([Inf], hospitalNumCap + 1),
    ones(hospitalNumCap + 1),
    Optim.Options(show_trace = false)
);

  5.900208 seconds (15.71 M allocations: 1.961 GiB, 6.81% gc time, 96.10% compilation time)

se = sqrt.(diag(inv(obj_for_Optim.H) ./ numObs));
DataFrame(
    estimates = optim_res.minimizer,
    se = se
)

7×2 DataFrame

Row	estimates	se
	Float64	Float64
1	54.6236	1.52851
2	33.5433	1.3618
3	8.08346	0.482284
4	4.23655	0.326509
5	2.51215	0.260241
6	1.63297	0.215186
7	1.39041	0.046226

5.3 推定値に基づくエクササイズ

alpha_est = optim_res.minimizer[1:hospitalNumCap];
gamma_est = optim_res.minimizer[hospitalNumCap + 1];

EntryThreshold = zeros(Int64, (hospitalNumCap, 2));

deno = alpha_est[1];
EntryThreshold[1,:] .= round(gamma_est / deno .* 1e+6);
for i = 2:hospitalNumCap
    deno = deno - alpha_est[i];
    EntryThreshold[i, :] = round.([gamma_est / deno .* 1e+6, gamma_est / deno .* 1e+6 / i]);
end

EntryThreshold    

6×2 Matrix{Int64}:
  25454  25454
  65958  32979
 106981  35660
 158717  39679
 222532  44506
 301269  50212