using CSV
using DataFrames
using StringEncodings
using FixedEffectModels
using RegressionTables
using Plots
using LinearAlgebra
using Statistics
using Optim
using Printf
using ForwardDiff
using Random
using GLM
using Serialization3 基礎編 2
mutable struct datalist_struct
X1::Array{Float64,2};
X2::Array{Float64,2};
Z::Array{Float64,2};
ShareVec::Vector{Float64};
marketIndex::Vector{Int64};
logitshare::Vector{Float64};
randomDrawMat::Array{Float64,2};
marketIndexMat::BitMatrix
end3.1 関数の読み込み
for file in readdir("functions/demand_estimation")
include("functions/demand_estimation/" * file)
end3.2 前回からの加工済みデータの読み込み
data = CSV.read("tmp/demand_estimation_1/data.csv", DataFrame);3.3 データの準備
NIPPYOautoIDvec = [
260, 4, 76, 104, 64, 54, 152, 153, 71, 197,
42, 45, 114, 208, 209, 77, 236, 58, 127, 187,
79, 175, 19, 117, 216, 112, 256, 119, 37, 158
];
dataNIPPYO = data[
in(NIPPYOautoIDvec).(data[:, :NameID]),
[:year, :share, :NameID, :Sales, :price, :hppw, :FuelEfficiency, :size, :Name]
];
dataNIPPYO[!, :log_sales] = log.(dataNIPPYO[:, :Sales]);
dataNIPPYO[!, :log_price] = log.(dataNIPPYO[:, :price]);
dataNIPPYO[!, :log10_sales] = log10.(dataNIPPYO[:, :Sales]);
dataNIPPYO[!, :log10_price] = log10.(dataNIPPYO[:, :price]);3.4 Logitモデルにおける価格弾力性行列の作成
data[!, :logit_share] = log.(data[:, :share]) .- log.(data[:, :share0]);
resultGH = reg(
data,
@formula(logit_share ~ (
price ~ iv_GH_own_hppw + iv_GH_own_FuelEfficiency + iv_GH_own_size +
iv_GH_other_hppw + iv_GH_other_FuelEfficiency + iv_GH_other_size
) + hppw + FuelEfficiency + size),
Vcov.robust(),
save = true
);
data2016 = dataNIPPYO[
dataNIPPYO.year .== 2016,
[:price, :share, :NameID, :Name]
];
numCars2016 = nrow(data2016);
ownElasticity = (
resultGH.coef[resultGH.coefnames .== "price"][1] .*
data2016.price .*
(1.0 .- data2016.share)
);
crossElasticity = (
- resultGH.coef[resultGH.coefnames .== "price"][1] .*
data2016.price .*
data2016.share
);
elasticityMat = reduce(
hcat,
[crossElasticity for car = 1:numCars2016]
);
elasticityMat[diagind(elasticityMat)] = ownElasticity;
elasticityMat[[12, 13, 10, 1], [12, 13, 10, 1]]4×4 Matrix{Float64}:
-1.76455 0.00114929 0.00114929 0.00114929
0.00122042 -0.818689 0.00122042 0.00122042
0.000148175 0.000148175 -0.959449 0.000148175
0.00184509 0.00184509 0.00184509 -0.671753.5 入れ子型ロジット推定のためのBLP操作変数の作成
transform!(
groupby(data, [:year, :Maker, :Type]),
(
[:hppw, :FuelEfficiency, :size] .=>
sum .=>
[:hppw_sum_own, :FuelEfficiency_sum_own, :size_sum_own]
),
(
[:hppw, :FuelEfficiency, :size] .=>
(x -> sum(x.^2)) .=>
[:hppw_sqr_sum_own, :FuelEfficiency_sqr_sum_own, :size_sqr_sum_own]
),
nrow => :group_n
);
transform!(
groupby(data, [:year, :Type]),
(
[:hppw, :FuelEfficiency, :size] .=>
sum .=>
[:hppw_sum_mkt, :FuelEfficiency_sum_mkt, :size_sum_mkt]
),
(
[:hppw, :FuelEfficiency, :size] .=>
(x -> sum(x.^2)) .=>
[:hppw_sqr_sum_mkt, :FuelEfficiency_sqr_sum_mkt, :size_sqr_sum_mkt]
),
nrow => :mkt_n
);
for variable in ["hppw", "FuelEfficiency", "size"]
addIVColumnsForNestedLogit!(data, variable)
end
data[!, :iv_BLP_own_num_nest] = data[:, :group_n] .- 1;
data[!, :iv_BLP_other_num_nest] = data[:, :mkt_n] .- data[:, :group_n];3.6 入れ子型ロジットモデルの推定
data = transform(
groupby(data, [:year, :Type]),
:Sales => sum => :sum_year_body
);
data[!, :inside_share] = data.Sales ./ data.sum_year_body;
data[!, :log_inside_share] = log.(data.Sales ./ data.sum_year_body);
resultOLS = reg(
data,
@formula(
logit_share ~
price + log_inside_share + hppw + FuelEfficiency + size
)
);
resultBLPNested = reg(
data,
@formula(logit_share ~ (
price + log_inside_share ~ iv_BLP_own_hppw_nest + iv_BLP_own_FuelEfficiency_nest + iv_BLP_own_size_nest +
iv_BLP_other_hppw_nest + iv_BLP_other_FuelEfficiency_nest + iv_BLP_other_size_nest +
iv_BLP_own_num_nest + iv_BLP_other_num_nest
) + hppw + FuelEfficiency + size),
Vcov.robust()
);
regtable(resultOLS, resultBLPNested)
-----------------------------------------------
logit_share
---------------------
(1) (2)
-----------------------------------------------
(Intercept) -7.557*** -9.548***
(0.166) (0.239)
price -0.307*** -0.654***
(0.013) (0.053)
log_inside_share 0.782*** 0.595***
(0.009) (0.035)
hppw 10.636*** 18.925***
(0.657) (1.965)
FuelEfficiency 0.055*** 0.069***
(0.004) (0.006)
size 0.156*** 0.227***
(0.008) (0.012)
-----------------------------------------------
Estimator OLS IV
-----------------------------------------------
N 1,823 1,823
R2 0.862 0.765
Within-R2
First-stage F statistic 1.667
-----------------------------------------------
3.6.1 入れ子型ロジットモデルにおける自己価格弾力性の計算
alpha1 = resultOLS.coef[resultOLS.coefnames .== "price"][1]
sigma1 = resultOLS.coef[resultOLS.coefnames .== "log_inside_share"][1]
alpha2 = resultBLPNested.coef[resultBLPNested.coefnames .== "price"][1]
sigma2 = resultBLPNested.coef[resultBLPNested.coefnames .== "log_inside_share"][1]
data[!, :own_elas_ols] = alpha1 .* data[:, :price] .* (
1.0 .- sigma1 .* data[:, :inside_share] .-
(1.0 .- sigma1) .* data[:, :share]
) ./ (1.0 .- sigma1);
data[!, :own_elas_ivblp_nested] = alpha2 .* data[:, :price] .* (
1.0 .- sigma2 .* data[:, :inside_share] .-
(1.0 .- sigma2) .* data[:, :share]
) ./ (1.0 .- sigma2);
describe(data[:, r"^own_elas"], :mean, :std, :median, :min, :max)2×6 DataFrame
| Row | variable | mean | std | median | min | max |
|---|---|---|---|---|---|---|
| Symbol | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | own_elas_ols | -3.52131 | 2.55849 | -2.84383 | -17.7936 | -0.96257 |
| 2 | own_elas_ivblp_nested | -4.04935 | 2.93775 | -3.2664 | -20.4015 | -1.11262 |
dataNIPPYO = data[
in(NIPPYOautoIDvec).(data[:, :NameID]),
[
:year, :share, :Type, :inside_share,
:NameID, :Sales, :price, :hppw, :FuelEfficiency, :size, :Name
]
];
dataNIPPYO[!, :log_sales] = log.(dataNIPPYO[:, :Sales]);
dataNIPPYO[!, :log_price] = log.(dataNIPPYO[:, :price]);
data2016 = dataNIPPYO[
dataNIPPYO.year .== 2016,
[:price, :Type, :share, :inside_share, :NameID, :Name]
];
numCars2016 = nrow(data2016);
ownElasticityNestedLogit = (
alpha2 .*
data2016.price .*
(
1.0 .- sigma2 .*
data2016.inside_share .-
(1.0 .- sigma2) .*
data2016.share
) ./
(1.0 .- sigma2)
);
crossElasticityOtherGroup = (
- alpha2 .*
data2016.price .*
data2016.share
);
crossElasticityOtherGroup = reduce(
hcat,
[crossElasticityOtherGroup for car = 1:numCars2016]
);
crossElasticityOtherGroup[
diagind(crossElasticityOtherGroup)
] = ownElasticityNestedLogit;crossElasticitySameGroup = reduce(
hcat,
[
- alpha2 .* data2016.price .* (
sigma2 .* data2016.inside_share .+ (1.0 .- sigma2) .* data2016.share
) ./ (1.0 .- sigma2)
for car in 1:numCars2016]
);
sameGroupIndicatorMat = (
reduce(hcat, [data2016.Type for car = 1:numCars2016]) .==
permutedims(reduce(hcat, [data2016.Type for car = 1:numCars2016]))
);
otherGroupIndicatorMat = (sameGroupIndicatorMat .== 0);
elasticityNestedLogitMat = (
crossElasticitySameGroup .*
sameGroupIndicatorMat .+
crossElasticityOtherGroup .*
otherGroupIndicatorMat
);
elasticityNestedLogitMat[
diagind(elasticityNestedLogitMat)
] = ownElasticityNestedLogit;
elasticityNestedLogitMat[[12, 13, 10, 1], [12, 13, 10, 1]]4×4 Matrix{Float64}:
-5.1197 0.0477573 0.0477573 0.00136172
0.050713 -2.34881 0.050713 0.001446
0.00615725 0.00615725 -2.80217 0.000175564
0.00218614 0.00218614 0.00218614 -1.832963.7 ランダム係数ロジットモデルの推定
3.7.1 下準備
sort!(data, [:year, :NameID]);
N = nrow(data);
T = length(unique(data.year));
X1 = hcat(ones(N), Matrix(data[:, [:price, :FuelEfficiency, :hppw, :size]]));
X2 = hcat(data.price, ones(N), data.size);
Z = hcat(
ones(N),
Matrix(data[:, [:FuelEfficiency, :hppw, :size]]),
Matrix(data[:, r"^iv_GH.*(?<!nest)$"])
);
Random.seed!(123);
Nsim = 500;randomDrawMat = randn(size(X2)[2], Nsim);
marketIndex = data.year;
uniqueMarketIndex = sort(unique(data.year));
marketIndexMat = reduce(
hcat,
[marketIndex .== market for market in uniqueMarketIndex]
);
datalist = datalist_struct(
X1,
X2,
Z,
data.share,
marketIndex,
data.logit_share,
randomDrawMat,
marketIndexMat
);3.7.2 推定
ここでの推定結果が元のサポートサイトのものとかなり異なることに注意されたい。 元の推定では、ランダム係数の推計値が初期値からほとんど変わっていないように見える。 もしかすると、本のコードでは市場シェアの更新（永続代入演算子が用いられている）が正しくされていないのかもしれない。
initial_x = [0.3, 18, 0.01];
delta_ini = calculateMeanUtil(initial_x, datalist, datalist.logitshare);objFunc_for_Optim = OnceDifferentiable(
x -> calculateGMMObjective(x, datalist, delta_ini),
initial_x;
autodiff = :forward
);
@time resultGMM = optimize(
objFunc_for_Optim,
[0.0, 0.0, 0.0],
[Inf, Inf, Inf],
initial_x,
Fminbox(LBFGS()),
Optim.Options(show_trace = false)
);resultGMM.minimizer3-element Vector{Float64}:
0.6407730207977352
14.72488892325358
1.689052459807443e-11W = inv(datalist.Z' * datalist.Z);
delta = calculateMeanUtil(resultGMM.minimizer, datalist, delta_ini);
beta_hat = (
(datalist.X1' * datalist.Z * W * datalist.Z' * datalist.X1) \
(datalist.X1' * datalist.Z * W * datalist.Z' * delta)
);
beta_hat5-element Vector{Float64}:
-31.984921285756943
-1.0841982623293536
0.1076677873797389
8.581036411670775
0.282984400128923233.7.3 標準誤差の計算
Xi = delta - X1 * beta_hat;
Omega_hat = reduce(+, Z[i,:] * Z[i,:]' .* Xi[i]^2 ./ N for i = 1:N);
Ddelta = ForwardDiff.jacobian(
x -> calculateMeanUtil(x, datalist, delta),
resultGMM.minimizer
);
G = Z' * hcat(- X1, Ddelta) ./ N;
AsyVarMat = (G' * W * G) \ G' * W * Omega_hat * W * G * inv(G' * W * G);
Ase = sqrt.(diag(AsyVarMat) ./ N);
DataFrame(
Var = [
"Const", "Price", "Fuel Efficiency", "hppw",
"size", "random_price", "random_constant", "random_size"
],
Est = vcat(beta_hat, resultGMM.minimizer),
se = Ase
)8×3 DataFrame
| Row | Var | Est | se |
|---|---|---|---|
| String | Float64 | Float64 | |
| 1 | Const | -31.9849 | 6.91814 |
| 2 | Price | -1.0842 | 1.34743 |
| 3 | Fuel Efficiency | 0.107668 | 0.0115801 |
| 4 | hppw | 8.58104 | 3.13587 |
| 5 | size | 0.282984 | 0.127824 |
| 6 | random_price | 0.640773 | 0.845573 |
| 7 | random_constant | 14.7249 | 21.1974 |
| 8 | random_size | 1.68905e-11 | 1.79294 |
3.8 価格弾力性行列の計算
market2016Index = data.year .== 2016;
elasmat = calculateElasticity(
data[market2016Index, :].price,
X2[market2016Index, :],
beta_hat,
resultGMM.minimizer,
randomDrawMat,
delta[market2016Index, :]
)
elasmat[[59, 80, 102, 113], [59, 80, 102, 113]]4×4 Matrix{Float64}:
-2.07694 0.0043804 0.00405532 0.00447937
0.0204111 -1.92602 0.0185025 0.0248292
0.00413169 0.00404559 -2.14446 0.00376725
0.00295524 0.0035155 0.00243949 -1.512163.9 応用: プライシング
以下で求められる値はサポートサイトのものとかなり異なるが、これは上で述べた推定結果が原因だと思われる。
price_range = range(1.8, 4.0, step = 0.05);
ownpi_res = calculateRevenue.(
price_range,
Ref(data), Ref(datalist), Ref(delta), Ref(beta_hat), Ref(resultGMM.minimizer),
"ownpi"
);
plot(price_range, ownpi_res, legend = false)
xlabel!("Price")
ylabel!("Revenue")price_range = range(1.8, 4.0, step = 0.05);
totalpi_res = calculateRevenue.(
price_range,
Ref(data), Ref(datalist), Ref(delta), Ref(beta_hat), Ref(resultGMM.minimizer),
"totalpi"
);
plot(price_range, totalpi_res, legend = false)
xlabel!("Price")
ylabel!("Revenue")ownpi_optim_res = optimize(
x -> - calculateRevenue(
x[1], data, datalist, delta,
beta_hat, resultGMM.minimizer, "ownpi"
),
[3.0]
);
@printf("Revenue-maximizing price: %.3f \n", ownpi_optim_res.minimizer[1])
@printf("Maximized revenue : %.3f", -ownpi_optim_res.minimum)Revenue-maximizing price: 3.256
Maximized revenue : 44811.735totalpi_optim_res = optimize(
x -> - calculateRevenue(
x[1], data, datalist, delta,
beta_hat, resultGMM.minimizer, "totalpi"
),
[3.0]
);
@printf("Revenue-maximizing price: %.3f \n", totalpi_optim_res.minimizer[1])
@printf("Maximized revenue : %.3f", -totalpi_optim_res.minimum)Revenue-maximizing price: 3.508
Maximized revenue : 570181.486calculateRevenue(
ownpi_optim_res.minimizer[1],
data,
datalist,
delta,
beta_hat,
resultGMM.minimizer,
"totalpi"
)569858.7798782628