-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathstatapp.m
272 lines (229 loc) · 10.2 KB
/
statapp.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
% 8/25/2021 Shuowen Chen and Hiroaki Kaido
% Computes the test statistic in application
function [Test, g_delta, varI, h, Vstar] = statapp(delta_hat, data, x1, x2)
% inputs:
% delta_hat: first-step estimated delta
% data: market outcomes
% x1: covariates of player 1, same definition for x2 (n by d)
% output:
% Test: test statstic
n = size(data, 1); % sample size
[zdelta, zbeta] = compute_z_allapp([0, 0], delta_hat);
outcome = count(data, x1, x2);
Cbeta = zbeta * outcome / sqrt(n);
Cdelta = zdelta * outcome / sqrt(n);
% Compute the outer product approximation of information matrix and
% decomposed to four parts for orthogonalization
% I_beta2 is 2 by 2; I_betadelta is 2 by 6, and I_delta2 is 6 by 6.
% The following is intermediate for I_beta2. Dim is 4 by 64, where 64
% denotes all possible combinations of market outcome and covariates
betablock = [zbeta(1,:).*zbeta(1,:); zbeta(1,:).*zbeta(2,:);...
zbeta(2,:).*zbeta(1,:); zbeta(2,:).*zbeta(2,:)];
% The following block is intermediate for I_delta2. Dim is 36 by 64.
% The first dimension is 36 since each player has 3 nuisance parameters and
% we need to consider each pairwise combination
deltablock = zeros(36, 64);
for i = 1:6
for j = 1:6
deltablock((j+6*(i-1)),:) = zdelta(i,:).*zdelta(j,:);
end
end
% The following block is intermediate for I_deltabeta and I_betadelta (which
% is the transpose of the former).
crossblock = zeros(12, 64);
for i = 1:2
for j = 1:6
crossblock((j+6*(i-1)),:) = zbeta(i,:).*zdelta(j,:);
end
end
% Muliply by outcome of each combination
betablock = betablock * outcome / n;
deltablock = deltablock * outcome / n;
crossblock = crossblock * outcome / n;
% Reshape
I_beta2 = reshape(betablock, 2, 2)'; % 2 by 2
I_delta2 = reshape(deltablock, 6, 6)'; % 6 by 6
I_betadelta = reshape(crossblock, 6, 2)'; % 2 by 6
% Regularization (barely matters for this application)
I_delta2 = get_sigmatilde(I_delta2);
g_delta = Cbeta - I_betadelta/I_delta2 * Cdelta;
varI = I_beta2 - I_betadelta/I_delta2 * I_betadelta';
varI = get_sigmatilde(varI);
% Compute the test statistic
gtilde = varI^(-0.5)*g_delta;
Test = gtilde'*gtilde;
Vmin = @(x) quadform(varI, g_delta-x);
[h, Vstar] = fmincon(Vmin, [-0.5; -0.5], eye(2), zeros(2, 1)); % solve the quadratic equation in the statistic
Test = Test - Vstar;
end
% Auxiliary functions
% 1. Count the frequencies of combination of outcome and covariates
function output = count(data, x1, x2)
% turn into binary code
% note we exclude the constants here since they always show up
temp = 10000*data + 1000*x1(:,1) + 100*x1(:,2) + 10*x2(:,1) + x2(:,2);
% turn into decimal number (0-63) so as to do counting using for loop
temp = bin2dec(num2str(temp));
% placeholder
output = NaN(64, 1);
for i = 0 : 63
output(i+1) = sum(temp == i);
end
end
% 2. Compute scores for each combination
function [z_delta, z_beta] = compute_z_allapp(beta, delta)
% Inputs:
% beta: structural parameter of interest (1 by 2 vector)
% delta: nuisance parameter (2*3 by 1)
% delta order: presence, size and constant
%
% Outputs:
% z_delta: a 6 by 64 matrix with the following form
% outcome (0,0) (0,1) (1,0) (1,1)
% covariate (0,0,1)(0,0,1) (0,0,1)(0,1,1) ...(1,1,1)(1,1,1)| | | |
% delta1_1 | | | |
% delta1_2 | | | |
% delta1_3 | | | |
% delta2_1 | | | |
% delta2_2 | | | |
% delta2_3 | | | |
%
% z_beta: a 2 by 64 matrix with the following form
% outcome (0,0) (0,1) (1,0) (1,1)
% covariate (0,0,1)(0,0,1) (0,0,1)(0,1,1) ...(1,1,1)(1,1,1)| | | |
% beta1_1 | | | |
% beta1_2 | | | |
%
% Placeholders:
z_delta = zeros(6, 64);
z_beta = zeros(2, 64);
% Construct the 4 possible covariate configurations
% Note: each player has the same 4 possible configurations, but when
% interacted there will be a totality of 16 cases
x = [0,0,1;0,1,1;1,0,1;1,1,1]; % constant is the 3rd column
% 4 possible combinations for x*delta for two players (4 by 2)
xdelta = [x*delta(1:3), x*delta(4:6)];
phi = normpdf(xdelta);
Phi = normcdf(xdelta);
phi_beta = normpdf(xdelta + beta);
Phi_beta = normcdf(xdelta + beta);
% For events (0, 0) and (1, 1) no subcase consideration
% zbeta_11 is 4 by 2
zbeta_11 = phi_beta./Phi_beta;
zbeta_00 = zeros(2, 16);
% Scores related to delta now features 4 possible combinations of
% covariates and 3 coefficients (each is 4 by 3)
zdelta1_00 = -phi(:,1).*x./(1-Phi(:,1));
zdelta2_00 = -phi(:,2).*x./(1-Phi(:,2));
zdelta1_11 = phi_beta(:,1).*x./Phi_beta(:,1);
zdelta2_11 = phi_beta(:,2).*x./Phi_beta(:,2);
% Need to change Phi, phi and x to accommodate the combined considerations
Phi1 = repelem(Phi(:,1),4);
Phi2 = repmat(Phi(:,2),[4,1]);
phi1 = repelem(phi(:,1),4);
phi2 = repmat(phi(:,2),[4,1]);
Phi_beta1 = repelem(Phi_beta(:,1),4);
Phi_beta2 = repmat(Phi_beta(:,2),[4,1]);
phi_beta1 = repelem(phi_beta(:,1),4);
phi_beta2 = repmat(phi_beta(:,2),[4,1]);
x_first = repelem(x,4,1);
x_second = repmat(x,[4,1]);
% For events (1, 0) and (0, 1), for each configuration of Xdelta,
% there are three subcases to consider. We have six conditions
% (two of which are repetitive) to consider the regions.
% Given each local alternative, we can evaluate the conditions and thus
% determine which subregion the given local alternative belongs to.
% First define the 5 variables used in the conditions
% Each is 16 by 1
var1 = Phi1.*(1-Phi_beta2);
z1 = Phi1.*(1-Phi2);
z2 = (1-Phi1).*Phi2 + Phi1 + Phi2 - Phi1.*Phi2 - Phi_beta1.*Phi_beta2;
var2 = (z1.*z2 - Phi2.*(1-Phi1).*z1)./(Phi1 + Phi2 - 2*Phi1.*Phi2);
var3 = Phi1.*(1-Phi2);
var4 = Phi_beta1.*Phi2;
var5 = Phi_beta1.*Phi_beta2;
% Pre-allocation. Note: each player's zdelta has 3 components
zdelta_0110 = zeros(6, 32);
zbeta_0110 = zeros(2, 32);
for i = 1:16 % loop over combined covariate configurations
% (0,1) chosen
if var4(i) + var3(i) - var5(i) > var2(i) && var1(i) > var3(i) + var4(i) - var5(i)
q1_10 = Phi1(i)*(1-Phi2(i))+Phi_beta1(i)*Phi2(i)-Phi_beta1(i)*Phi_beta2(i);
zbeta1_10 = (Phi2(i)*phi_beta1(i)-Phi_beta2(i)*phi_beta1(i))/q1_10;
zbeta2_10 = -Phi_beta1(i)*phi_beta2(i)/q1_10;
zdelta1_10 = x_first(i,:) * (phi1(i)*(1-Phi2(i))+Phi2(i)*phi_beta1(i)-...
Phi_beta2(i)*phi_beta1(i))/q1_10;
zdelta2_10 = x_second(i,:) * (-phi2(i)*Phi_beta1(i)+Phi1(i)*phi2(i)-...
Phi_beta1(i)*phi_beta2(i))/q1_10;
zdelta1_01 = -x_first(i,:)*phi_beta1(i)/(1-Phi_beta1(i));
zdelta2_01 = x_second(i,:)*phi2(i)/Phi2(i);
zbeta1_01 = -phi_beta1(i)/(1-Phi_beta1(i));
zbeta2_01 = 0;
% (1,0) chosen
elseif var1(i) < var2(i) && var1(i) - var3(i) - var4(i) + var5(i) > 0
q1_01 = (1-Phi1(i))*Phi2(i) + Phi_beta2(i)*(Phi1(i)-Phi_beta1(i));
zdelta1_10 = x_first(i,:)*phi1(i)/Phi1(i);
zdelta2_10 = -x_second(i,:)*phi_beta2(i)/(1-Phi_beta2(i));
zbeta1_10 = 0;
zbeta2_10 = -phi_beta2(i)/(1-Phi_beta2(i));
zdelta1_01 = (-x_first(i,:)*Phi2(i)*phi1(i) + Phi_beta2(i)*...
x_first(i,:)*(phi1(i)-phi_beta1(i)))/q1_01;
zdelta2_01 = (x_second(i,:)*(1-Phi1(i))*phi2(i) + x_second(i,:)*...
(Phi1(i)-Phi_beta1(i))*phi_beta2(i))/q1_01;
zbeta1_01 = -Phi_beta2(i)*phi_beta1(i)/q1_01;
zbeta2_01 = (Phi1(i)-Phi_beta1(i))*phi_beta2(i)/q1_01;
else
% Mixture over (1,0) and (0,1)
denominator1 = Phi1(i)+Phi2(i)-Phi1(i)*Phi2(i)-Phi_beta1(i)*Phi_beta2(i);
denominator2 = Phi1(i)+Phi2(i)-2*Phi1(i)*Phi2(i);
zdelta1_10 = phi1(i)*x_first(i,:)/Phi1(i) + (phi1(i)*x_first(i,:)...
-phi1(i)*x_first(i,:)*Phi2(i)-phi_beta1(i)*Phi_beta2(i)*...
x_first(i,:))/denominator1 - (phi1(i)*x_first(i,:)*(1-2*Phi2(i)))/denominator2;
zdelta2_10 = -phi2(i)*x_second(i,:)/(1-Phi2(i)) + (phi2(i)*...
x_second(i,:)-phi2(i)*Phi1(i)*x_second(i,:)-phi_beta2(i)*...
Phi_beta1(i)*x_second(i,:))/denominator1 - (phi2(i)*x_second(i,:)*...
(1-2*Phi1(i)))/denominator2;
zbeta1_10 = -phi_beta1(i)*Phi_beta2(i)/denominator1;
zbeta2_10 = -Phi_beta1(i)*phi_beta2(i)/denominator1;
zdelta1_01 = -phi1(i)*x_first(i,:)/(1-Phi1(i)) + (phi1(i)*...
x_first(i,:)-phi1(i)*Phi2(i)*x_first(i,:)-phi_beta1(i)*...
Phi_beta2(i)*x_first(i,:))/denominator1 - (phi1(i)*...
x_first(i,:)*(1-2*Phi2(i)))/denominator2;
zdelta2_01 = phi2(i)*x_second(i,:)/Phi2(i) + (phi2(i)*x_second(i,:)-...
phi2(i)*Phi1(i)*x_second(i,:)-phi_beta2(i)*Phi_beta1(i)*...
x_second(i,:))/denominator1 - (phi2(i)*x_second(i,:)*(1-2*Phi1(i)))/denominator2;
zbeta1_01 = -phi_beta1(i)*Phi_beta2(i)/denominator1;
zbeta2_01 = -phi_beta2(i)*Phi_beta1(i)/denominator1;
end
% stacking up the results
zbeta_0110(:, 16+i) = [zbeta1_01; zbeta2_01];
zbeta_0110(:, 32+i) = [zbeta1_10; zbeta2_10];
zdelta_0110(:, 16+i) = [zdelta1_01'; zdelta2_01'];
zdelta_0110(:, 32+i) = [zdelta1_10'; zdelta2_10'];
end
% Fill in the rest of z_beta
z_beta(:,1:16) = zbeta_00;
z_beta(:,49:end) = repelem(zbeta_11', 1, 4);
% Fill in the rest of z_delta
z_delta(1:3,1:16) = repelem(zdelta1_00', 1, 4);
z_delta(4:end,1:16) = repelem(zdelta2_00', 1, 4);
z_delta(1:3,49:end) = repelem(zdelta1_11', 1, 4);
z_delta(4:end,49:end) = repelem(zdelta2_11', 1, 4);
end
% 3. Regularization using Andrews and Barwick
function SigmaTilde = get_sigmatilde(Sigma)
% Inputs
% Sigma: Covariance matrix
% Output
% SigmaTilde: Regularized covariance matrix based on Andrews and Barwick
sh_param = 0.05;
sigma = sqrt(diag(Sigma));
Omega = Sigma./(sigma*sigma');
D = diag(diag(Sigma));
SigmaTilde = Sigma + max(sh_param-det(Omega), 0).*D;
end
% 4. Define quadratic form
function t = quadform(A, g)
gtilde = A^(-0.5)*g;
t = gtilde'*gtilde;
end