-
Notifications
You must be signed in to change notification settings - Fork 4
/
boottest.ado
1262 lines (1144 loc) · 57.6 KB
/
boottest.ado
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
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
*! boottest 4.4.11 12 April 2024
*! Copyright (C) 2015-24 David Roodman
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
* Version history at bottom
cap program drop boottest
program define boottest
version 11
cap version 13.1
if _rc {
di as err "This version of {cmd:boottest} requires Stata version 13 or later. An older version compatible with Stata `c(stata_version)'"
di as err "is at https://github.com/droodman/boottest/releases/tag/v2.6.0."
exit _rc
}
cap noi _boottest `0'
local rc = _rc
if `'"`env'"'!="" qui jl SetEnv `env' // revert to original Julia environment local env set by _boottest
constraint drop `anythingh0'
cap mata mata drop _boottestp
cap mata mata drop _boottestC
cap mata mata drop _boottestkeepC
cap mata mata drop _boottestm
cap mata mata drop _boottestt
cap drop `contourvars'
exit `rc'
end
cap program drop _boottest
program define _boottest, rclass sortpreserve
version 11
local JLVERSION 0.8.0
local cmd = cond(substr("`e(cmd)'", 1, 6)=="ivreg2" | ("`e(cmd)'"=="ivreghdfe" & "`e(extended_absvars)'"==""), "ivreg2", "`e(cmd)'")
local ivcmd = cond(inlist("`cmd'","reghdfe","ivreghdfe"), cond("`e(model)'"=="iv", "ivreg2", ""), cond("`cmd'"=="xtivreg2", "ivreg2", "`cmd'"))
if "`e(cmd)'" == "" {
di as err "No estimates detected."
error 198
}
if "`e(prefix)'" == "svy" {
di as err "Doesn't work after {cmd:svy}."
exit 198
}
if inlist("`cmd'", "mvreg", "sureg") {
di as err "Doesn't work after {cmd:`e(cmd)'}."
exit 198
}
if "`e(cmd)'" == "margins" {
di as err "Doesn't work after {cmd:margins ..., post}."
exit 198
}
if inlist("`cmd'", "xtreg", "xtivreg") & "`e(model)'"!="fe" {
di as err "Doesn't work after {`cmd', `e(model)'}."
exit 198
}
if "`cmd'"=="xtivreg2" & "`e(xtmodel)'"!="fe" {
di as err "Doesn't work after {`cmd', `e(xtmodel)'}."
exit 198
}
if inlist("`cmd'","reghdfe","ivreghdfe","reghdfejl","ivreghdfejl") & `"`e(absvars)'"'!="" {
fvunab absvars: `e(absvars)'
if `:word count `absvars''>1 {
di as err "Doesn't work after {cmd:`cmd'} with more than one set of absorbed fixed effects or with absorbed interaction terms."
exit 198
}
if strpos("`absvars'", "c.") {
di as err "Doesn't work after {cmd:`cmd'} with absorbed interaction terms containing slopes."
exit 198
}
else local absvars = subinstr(subinstr("`absvars'", "#", " ", .), "i.", "", .)
}
if inlist("`cmd'", "sem", "gsem") & c(stata_version) < 14 {
di as err "Requires Stata version 14.0 or later to work with {cmd:`e(cmd)'}."
exit 198
}
if "`e(cmd)'"=="regress" & "`e(model)'" == "iv" {
di as err "Doesn't support the undocumented IV syntax of {cmd:regress}."
exit 198
}
if ("`ivcmd'"=="ivreg2" & "`e(model)'"=="gmm2s") | ("`cmd'"=="ivregress" & "`e(estimator)'"=="gmm") {
di as err "GMM no longer supported."
exit 198
}
tokenize `"`0'"', parse(",")
if `"`1'"' != "," {
local h0s `1'
macro shift
}
local 0 `*'
syntax, [h0(numlist integer >0) Reps(integer 999) seed(string) BOOTtype(string) CLuster(string) Robust BOOTCLuster(string) noNULl QUIetly WEIGHTtype(string) Ptype(string) STATistic(string) NOCI Level(real `c(level)') NOSMall SMall SVMat ///
noGRaph gridmin(string) gridmax(string) gridpoints(string) graphname(string asis) graphopt(string asis) ar MADJust(string) CMDline(string) MATSIZEgb(real 1000000) PTOLerance(real 1e-3) svv MARGins ///
issorted julia PRECision(integer 64) Format(string) jk JACKknife *]
if "`format'"=="" local format %10.4g
qui query born
local mlabformat = cond($S_1 < td(16oct2019), "", `"mlabformat(`format')"')
local jk = "`jk'`jackknife'" != ""
local quietly = "`quietly'" != ""
local julia = "`julia'" != ""
if `julia' & !0$boottest_julia_loaded {
cap jl version
if _rc {
di as err `"Can't access Julia. {cmd:boottest} requires that the {cmd:jl} command be installed, via {stata ssc install julia}."
di as err "And it requires that Julia be installed, following the instruction under Installation in {help jl##installation:help jl}."
exit 198
}
parse "`r(version)'", parse(".")
local v1: copy local 1
local v2: copy local 3
local v3: copy local 5
parse "`JLVERSION'", parse(".")
local version: di %05.0f `v1' %05.0f `v2' %05.0f `v3'
if "`version'" < "000010000000000" {
di as txt "The Stata package {cmd:julia} is not up to date. Attempting to update it with {stata ssc install julia, replace}." _n
ssc install julia, replace
}
qui jl GetEnv
c_local env `r(env)'
qui jl SetEnv boottest
jl AddPkg StableRNGs
jl AddPkg WildBootTests, minver(0.9.12)
_jl: using StableRNGs, WildBootTests;
_jl: rng = StableRNG(0); // create now, seed later
global boottest_julia_loaded 1
}
if "`small'" != "" & "`nosmall'" != "" {
di as err "{cmd:small} and {cmd:nosmall} options conflict."
exit 198
}
local margins = "`margins'" != ""
if `margins' & `"`h0s'`h0'"' != "" {
di as err "Include the {cmd:margins} option or state null hyptheses, but don't do both."
exit 198
}
if `margins' & "`r(predict1_label)'" != "Linear prediction" {
di as err _n "{cmd:margins} option only works with linear margins predictions such as average predictions after {cmd:regress}."
exit 198
}
if inlist("`e(cmd)'", "didregress", "xtdidregress") & `"`h0s'`h0'"' == "" local h0s r1vs0.`e(treatment)'
if `matsizegb'==1000000 local matsizegb .
if "`svv'" != "" tempname svv
if `reps'==0 local svmat
else {
local _svmat = cond("`svmat'"=="", "", "t")
local 0, `options'
syntax, [SVMat(string) *]
if !inlist(`"`svmat'"', "numer", "t", "") {
di as err "option " as res "svmat(" `svmat' ")" as err " not allowed."
error 198
}
if "`svmat'" == "" local svmat: copy local _svmat
}
if `"`options'"' != "" {
if `"`options'"'=="ci" di as err "Option {cmd:ci} is obsolete, because it is now the default."
else di as err `"Option `options' not allowed."'
exit 198
}
if `reps' < 0 {
di as err "{cmdab:r:eps()} option must be non-negative."
exit 198
}
local 0, `madjust'
syntax, [Bonferroni Sidak]
local madjust `bonferroni'`sidak'
if `"`graphname'"' != "" {
local 0: copy local graphname
syntax [anything(name=graphname)], [replace]
}
else local graphname Graph
if `"`gridpoints'"'!="" {
foreach g of local gridpoints {
cap confirm integer number `g'
local rc = _rc
if !_rc {
local rc = `g' <= 0
}
if `rc' {
di as err "{cmd:gridpoints()} entry not a positive integer."
exit 198
}
}
}
foreach macro in gridmin gridmax {
if `"``macro''"'!="" {
foreach g of local `macro' {
if `"`g'"' != "." {
cap confirm number `g'
if _rc {
di as err "{cmd:`macro'()} entry not numeric."
exit 198
}
}
}
}
}
if `"`robust'`cluster'"' != "" {
local hasrobust 1
local clustvars: copy local cluster
if `"`clustvars'"'!="" confirm var `clustvars'
local override 1
di as txt _n "Overriding estimator's cluster/robust settings with " as res `"`=cond("`clustvars'"=="", "robust", "cluster(`clustvars')")'"'
}
else {
local clustvars `e(clustvar)'
local hasrobust = "`e(vcetype)'" == "Robust" | "`clustvars'" != ""
}
local hasclust = `"`clustvars'"' != ""
local wtype `e(wtype)'
local null = "`null'" == ""
if "`svmat'"!="" tempname dist
local ar = "`ar'" != ""
if `ar' & `"`h0s'`h0'"' == "" local h0s `e(instd)'
if `:word count `weighttype'' > 1 {
di as err "The {cmd:weight:type} option must be {cmdab:rad:emacher}, {cmdab:mam:men}, {cmdab:nor:mal}, {cmdab:web:b}, or {cmdab:gam:ma}."
exit 198
}
if "`weighttype'"'=="" local weighttype rademacher
else {
local 0, `weighttype'
syntax, [RADemacher MAMmen NORmal WEBb GAMma]
local weighttype `rademacher'`mammen'`normal'`webb'`gamma'
}
if `:word count `ptype'' > 1 {
di as err "The {cmd:wp:type} option must be {cmdab:sym:metric}, {cmdab:eq:qualtail}, {cmd:lower}, or {cmd:upper}."
exit 198
}
if "`ptype'"'=="" local ptype symmetric
else {
local 0, `ptype'
syntax, [SYMmetric EQualtail LOWer UPper]
local ptype `symmetric'`equaltail'`lower'`upper'
}
if `"`statistic'"'=="" local statistic t
else if !inlist(`"`statistic'"', "t", "c") {
di as err "The {cmd:stat:istic} option must be {cmd:t} or {cmd:c}."
exit 198
}
else if "`statistic'"=="c" & `reps'==0 {
di as err "{cmdab:stat:istic(c)} not allowed with non-bootstrap tests."
exit 198
}
local ML = e(converged) != . & !inlist("`cmd'", "reghdfe", "ivreghdfe", "reghdfejl")
local IV = "`e(instd)'`e(endogvars)'" != ""
local LIML = ("`cmd'"=="ivreg2" & "`e(model)'"=="liml") | ("`cmd'"=="ivregress" & "`e(estimator)'"=="liml") | (inlist("`cmd'","reghdfe","ivreghdfe") & strpos("`e(title)'", "LIML"))
local WRE = `"`boottype'"'!="score" & `IV' & `reps'
local small = (e(df_r) != . | "`small'" != "" | e(cmd)=="cgmreg") & "`nosmall'"==""
local partial = `:word count `e(partial1)'' & inlist("`e(cmd)'", "ivreg2", "xtivreg2", "ivreghdfe")
local fuller `e(fuller)' // "" if missing
local K = e(kclass) // "." if missing
local DID = inlist("`cmd'", "didregress", "xtdidregress")
if `DID' {
if "`e(aggmethod)'" != "" {
di as err "Doesn't work after {cmd:`e(cmd)'} with aggregation."
exit 198
}
local treatment `e(treatment)'
local 0 `e(cmdline)'
syntax [anything], [NOGTEFFECTS *]
local DID = "`nogteffects'" == ""
}
local tz = cond(`small', "t", "z")
local symmetric Prob>|`tz'|
local equaltail 2 * min(Prob>|`tz'|, Prob<-|`tz'|)
local lower Prob<`tz'
local upper Prob>`tz'
if `ar' & !`IV' {
di as err "Anderson-Rubin test is only for IV models."
exit 198
}
if `jk' & `ML' {
di as err "boottest can't jackknife ML-based estimates."
exit 198
}
if "`boottype'"'=="" local boottype = cond(`ML', "score", "wild")
else {
local 0, `boottype'
syntax, [Wild SCore]
local boottype `score'`wild'
if "`boottype'" == "wild" & `ML' {
di as err "{cmd:boottype(wild)} not accepted after Maximum Likelihood-based estimation."
exit 198
}
if "`boottype'"=="score" & "`fuller'`K'" != "." & "`e(model)'"!="liml" {
di as err "{cmd:boottype(score)} not accepted after Fuller LIML or k-class estimation."
exit 198
}
}
local scoreBS = "`boottype'"=="score"
if inlist("`cmd'","xtreg","xtivreg","xtivreg2") local NFE = e(N_g) + 0`e(singleton)'
else local NFE = cond(`DID', 0`e(N_clust)', ///
cond("`cmd'"=="areg", 1+e(df_a), ///
max(0`e(K1)', 0`e(df_a)', 0`e(df_a_initial)'))) // reghdfe
local _FEname = cond(inlist("`cmd'","xtreg","xtivreg","xtivreg2"), "`e(ivar)'", cond(`DID', "`e(clustvar)'", "`e(absvar)'`absvars'"))
if `"`_FEname'"' != "" {
cap confirm numeric var `_FEname'
if _rc | `:word count `_FEname' > 1' {
tempvar FEname
qui egen long `FEname' = group(`_FEname') if e(sample)
}
else local FEname: copy local _FEname
}
if "`cmd'" == "xtreg" {
local 0 `e(cmdline)'
syntax [anything] [if] [in] [fw aw pw iw], [dfadj *]
local FEdfadj = ("`dfadj'" != "") * `NFE'
}
else if inlist("`cmd'","xtivreg","xtivreg2","xtdidregress") local FEdfadj 0
else if inlist("`cmd'", "reghdfe", "ivreghdfe") local FEdfadj = max(1, e(df_a))
else local FEdfadj: copy local NFE
if `"`seed'"'!="" set seed `seed'
tempname p padj se teststat df df_r hold C1 C R1 R r1 r1r R1R r b V b0 V0 keepC repsname repsFeasname t NBootClustname marginsH0 touse
mat `b' = e(b)
if "`e(wtype)'" != "" {
tokenize `e(wexp)', parse("=")
cap confirm var `2'
if _rc {
tempname wtname
qui gen double `wtname' `e(wexp)'
}
else if c(varabbrev)=="on" unab wtname: `2'
else local wtname: copy local 2
}
if `"`h0s'"' != "" {
if "`h0'" != "" {
di as err "Specify hypotheses before comma or with {cmd:h0()} option, but not both."
exit 198
}
local multiple = strpos(`"`h0s'"', "{")
if !`multiple' local h0s {`h0s'}
local N_h0s 0 // number of nulls
while `"`h0s'"' != "" {
gettoken h0 h0s: h0s, parse("{}")
if !inlist(`"`h0'"', "{", "}") {
local ++N_h0s
while `"`h0'"' != "" {
gettoken cns h0: h0, parse("()") match(m)
constraint free
constraint `r(free)' `cns'
local h0_`N_h0s' `h0_`N_h0s'' `r(free)'
local constraints `constraints' `r(free)'
c_local anythingh0 `constraints'
}
}
}
local anythingh0 1
}
else if `margins' {
if strlen("`r(Jacobian)'") {
mata _boottestp = selectindex(rowsum(st_matrix("r(Jacobian)"):!=0)) // skip all-zero rows
mata st_local("marginsnames", invtokens(st_matrixrowstripe("r(Jacobian)")[_boottestp,2]'))
mata st_matrix("`marginsH0'", st_matrix("r(Jacobian)")[_boottestp,])
mata st_local("N_h0s", strofreal(length(_boottestp)))
if `N_h0s'==0 {
di as err "No valid {cmd:margins} results not found."
error 303
}
scalar `df' = 1 // always when working with margins results, df = 1 and constant term = 0
mat `r' = 0
local 0 r(cmdline)
marksample marginstouse, strok
}
else {
di as err "{cmd:margins} results not found."
error 303
}
}
else {
local N_h0s 1 // number of nulls
if "`h0'" == "" {
di as txt _n "({cmd:h0(1)} assumed)"
local h0 1
}
if !`quietly' {
foreach c of numlist `h0' {
if `"`: constraint `c''"' == "" di as res "Constraint `c' not found and will be skipped."
}
}
local h0_1: copy local h0
}
if `N_h0s'==1 local madjust
makecns
if "`e(Cns)'" != "" {
local hascns 1
mat `C1' = e(Cns)
mat `R1' = `C1'[1...,1..colsof(`C1')-1]
mat `r1' = `C1'[1..., colsof(`C1') ]
}
cap _estimates drop `hold'
if !`ML' {
if ("`ivcmd'"=="ivreg2" & !inlist("`e(model)'", "ols", "iv", "gmm2s", "liml")) | ("`ivcmd'"=="ivregress" & !inlist("`e(estimator)'", "liml", "2sls", "gmm")) {
di as err "Only for use with OLS, 2SLS, single-equation GMM, and ML-based estimators."
exit 198
}
local Ynames `e(depvar)'
local colnames: colnames e(b)
if `partial' local colnames `colnames' `e(partial1)'
local k = `:word count `colnames'' + (`partial' & 0`e(partialcons)')
local _cons _cons
local cons = "`:list colnames & _cons'" != ""
local colnames: list colnames - _cons
if `IV' {
local Xnames_endog `e(instd)'
local Xnames_exog: list colnames - Xnames_endog
local cons = `cons' | e(partialcons)==1
if inlist("`cmd'", "ivreg2", "xtivreg2", "ivreghdfe") local ZExclnames `e(exexog1)'
// else if "`e(cmd)'"=="ivreghdfe" { // ivreghdfe provides no collinear instrument info
// if "`FEname'"=="" local ZExclnames `e(exexog1)'
// else {
// qui _rmdcoll `Ynames' `Xnames_exog' i.`FEname' `e(exexog1)' if e(sample) [`e(wtype)'`e(wexp)'] // will fail if number of groups exceeds matsize
// local varlist `r(varlist)'
// local n: word count `varlist'
// forvalues i=`=`n'-`:word count `e(exexog1)''+1'/`n' {
// local var: word `i' of `varlist'
// _ms_parse_parts `var'
// if !r(omit) local ZExclnames `ZExclnames' `var'
// }
// }
// }
else {
local ZExclnames `e(insts)'
local ZExclnames: list ZExclnames - Xnames_exog
}
}
else {
local Xnames_exog: copy local colnames
if inlist("`e(cmd)'", "didregress", "xtdidregress") local Xnames_exog: subinstr local Xnames_exog "r1vs0." ""
}
local cons = `cons' & !0`NFE' // no constant in FE model
if !`cons' local _cons
mata _boottestp = J(`cons', 1, `k') \ order(tokens("`colnames'")', 1)[invorder(order(tokens("`Xnames_exog' `Xnames_endog'")', 1))] // for putting vars in cons-exog-endog order
if `cons' mat `keepC' = 1
local colnames `_cons' `Xnames_exog' `Xnames_endog'
foreach varlist in Xnames_exog Ynames Xnames_endog ZExclnames {
fvrevar ``varlist'' if e(sample)
local revarlist `r(varlist)'
local _revarlist
forvalues i=1/`:word count ``varlist''' {
local var: word `i' of ``varlist''
_ms_parse_parts `var'
if !r(omit) {
local _revarlist `_revarlist' `: word `i' of `revarlist''
if inlist("`varlist'", "Xnames_exog", "Xnames_endog") {
local pos: list posof "`var'" in colnames
if `pos' mat `keepC' = nullmat(`keepC'), `pos'
}
}
}
local `varlist': copy local _revarlist
}
mata _boottestkeepC = st_matrix("`keepC'"); _boottestp = cols(_boottestkeepC)? _boottestp[_boottestkeepC] : J(0,1,0)
if 0`hascns' {
if `partial' mat `R1' = `R1', J(rowsof(`R1'), `:word count `e(partial1)'', 0) // add blanks for partialled-out
mata st_matrix("`R1'", st_matrix("`R1'")[,_boottestp]) // put cols in standardized order & drop those for omitted regressors
mata _boottestkeepC = rowsum(st_matrix("`R1'"):!=0)
mata st_matrix("`R1'", select(st_matrix("`R1'"), _boottestkeepC)) // drop rows corresponding purely to omitted variables
mata st_matrix("`r1'", select(st_matrix("`r1'"), _boottestkeepC))
}
if `cons' local Xnames_exog `hold' `Xnames_exog' // add constant term
}
local NErrClustVar : word count `clustvars'
if `hasclust' {
if 0`override' {
cap assert !missing(`:subinstr local clustvars " " ",", all') if e(sample), `=cond(c(stata_version)>=12.1, "fast", "")'
if _rc {
di as err "A clustering variable has a missing value for at least one observation."
exit 9
}
}
if `"`bootcluster'"' == "" {
local bootcluster: copy local clustvars
if `reps' & `NErrClustVar'>1 di as txt "({cmdab:bootcl:uster(`clustvars')} assumed)"
}
local clustvars `:list clustvars & bootcluster' `:list clustvars - bootcluster'
local allclustvars `:list bootcluster - clustvars' `clustvars' // put bootstrapping clusters first, error clusters last, overlap in middle
if "`issorted'"=="" sort `clustvars' `:list bootcluster - clustvars', stable // bootstrapping-only clusters left out of this sort. Effect is that in subcluster bootstrap, bootstrapping clusters are intersections of all clusterings
foreach clustvar in `allclustvars' {
if !inlist("`:type `clustvar''", "byte", "int", "long") {
tempvar IDname
qui egen long `IDname' = group(`clustvar') if e(sample)
local _allclustvars `_allclustvars' `IDname'
}
else local _allclustvars `_allclustvars' `clustvar'
}
local allclustvars: copy local _allclustvars
}
local NBootClustVar: word count `bootcluster'
forvalues h=1/`N_h0s' { // loop over multiple independent constraints
if `margins' mat `R' = `marginsH0'[`h', 1...]
else {
_estimates hold `hold', restore
ereturn post `b'
mat `b' = e(b)
local h0text
foreach c in `h0_`h'' {
local h0text = `"`h0text'"' + `" "`: constraint `c''""'
}
qui makecns `h0_`h'' // process hypothesis constraints into e(Cns)
if "`h0_`h''" != "`r(clist)'" {
local clist `r(clist)'
local clist: list h0_`h' - clist
foreach c in `clist' {
di as txt `"note: constraint `=cond(0`anythingh0', `"{res}`:constraint `c''{txt}"', "number `c'")' caused error {search r(111)}"'
if `DID' di `"If your test relates to the treatment effect, you may need to reference "{cmd:r1vs0.`treatment'}" instead of "{cmd:`treatment'}"."' _n
}
exit 111
}
cap mat `C' = e(Cns)
if _rc exit 111
_estimates unhold `hold'
if r(k_autoCns) mat `C' = `C'[r(k_autoCns)+1...,1...]
mat `R' = `C'[1...,1..colsof(`C')-1]
mat `r' = `C'[1..., colsof(`C') ]
* get rid of any remaining o. and b. constraints; r(k_autoCns) doesn't capture all
mata _boottestC = J(1,0,0)
foreach var in `:colnames `R'' {
_ms_parse_parts `var'
mata _boottestC = _boottestC, !`r(omit)'
}
mata _boottestC = rowsum(select(st_matrix("`R'"), _boottestC) :!= 0)
mata st_matrix("`R'", select(st_matrix("`R'"), _boottestC))
mata st_matrix("`r'", select(st_matrix("`r'"), _boottestC))
cap mat `R'[1,1] = `R'[1,1]
if _rc {
di as error _n "Null constraint applies only to omitted variables or base levels of factor variables."
continue
}
scalar `df' = rowsof(`R')
}
if 0`NFE' {
mata st_local("rc", strofreal(any(0:!=select(st_matrix("`R'"),st_matrixcolstripe("`R'")[,2]':=="_cons"))))
if `rc' {
di as err "In fixed-effect models, null hypotheses may not involve constant term."
exit 111
}
}
local plotmat
local peakmat
local cimat
if "`noci'"=="" & !`ML' & (`df'<=1 | `df'==2 & "`graph'"=="") {
if `reps' | !`scoreBS' | `null' | "`graph'"=="" tempname plotmat peakmat // don't request plot if doing simple Wald with no graph
if `level'<100 & `df'==1 tempname cimat
}
if `df'>1 & "`ptype'"!="symmetric" di as txt "Note: {cmd:ptype(`ptype')} ignored for multi-constraint null hypotheses."
if `df'<=2 { // a bit duplicative of code in the if-`ML' block just below...
if "`graph'"=="" { // construct axis label(s)
if `margins' local constraintLHS1 "Margins of `:word `h' of `marginsnames''"
else {
local coleq: coleq `b'
if "`:word 1 of `coleq''"=="_" local coleq
local colnames: colnames `b'
forvalues i=1/`=`df'' {
local terms 0
local constraintLHS`i'
forvalues j=1/`=colsof(`R')' {
if `R'[`i',`j'] {
if `terms++' local constraintLHS`i' `constraintLHS`i''+
local constraintLHS`i' `constraintLHS`i'' `=cond(`R'[`i',`j']!=1,"`=`R'[`i',`j']'*","")'`=cond("`coleq'"=="","","[`:word `j' of `coleq'']")'`:word `j' of `colnames''
}
}
}
}
}
foreach option in gridmin gridmax gridpoints {
if "``option''" == "" local `option' = `df' * ". "
else if `df' != `:word count ``option''' {
di as err "{cmd:`option'()} option needs " cond(`df'==1, "one entry", "two entries")
exit 199
}
tempname `option'vec
mata st_matrix("``option'vec'", strtoreal(tokens("``option''")))
}
}
if `ML' {
local K .
local k = colsof(`b')
if `null' {
if "`e(cmdline)'"=="" & `"`cmdline'"'=="" {
di as err "Original estimation command line needed. Include it in a {cmd:cmdline()} option."
exit 198
}
if "`cmd'"=="tobit" & c(stata_version)<15 {
di as err "Cannot impose null after {help tobit} in Stata versions below 15. Try re-fitting the model with {stata ssc describe cmp:cmp}, {help intreg}, or {help gsem}."
exit 198
}
mat `R1R' = `R' \ nullmat(`R1') // add null to model constraints
mat `r1r' = `r' \ nullmat(`r1')
if !`quietly' di as res _n "Re-running regression with null imposed." _n
if `"`cmdline'"'=="" local cmdline `e(cmdline)'
_parse expand lc lg: cmdline, common(CONSTraints(passthru) from(passthru) INIt(passthru) ITERate(passthru) Robust) // would be nice to extract, thus remove, cluster() option for speed, but it affects sample
local 0, `lg_op'
syntax, [from(passthru) INIt(passthru) ITERate(passthru) *]
local 0: copy local lc_1
syntax [anything] [aw pw fw iw] [if] [in], [*]
if "`e(cmd)'"=="ml" local max max // tack on max option in case ml run in interactive model and cmdline() is missing it
cap _estimates drop `hold'
_estimates hold `hold', restore
local coleq: coleq `b'
if "`:word 1 of `coleq''"=="_" local coleq
local colnames: colnames `b'
local _constraints
forvalues i=1/`=rowsof(`R1R')' {
local terms 0
local _constraint
forvalues j=1/`=colsof(`R1R')' {
if `R1R'[`i',`j'] {
if `terms++' local _constraint `_constraint' +
local eq = cond("`coleq'"=="" | inlist("`cmd'","xttobit","xtintreg") & "`coleq'"!=".", "", "[`:word `j' of `coleq'']") // hack around xttobit/xtintreg bug internally renaming first eq to "eq1"
local _constraint `_constraint' `=cond(`R1R'[`i',`j']!=1,"`=`R1R'[`i',`j']' * ","")' `eq'`:word `j' of `colnames''
}
}
local _constraint `_constraint' = `=`r1r'[`i',1]'
constraint free
local _constraints `_constraints' `r(free)'
constraint `r(free)' `_constraint'
}
local cmdline version `c(stata_version)': `anything' `if' `in' `=cond("`weight'"=="","",`"[`weight'`exp']"')', `max' `options' // "if `hold'" would prematurely restrict sample, interact badly with probit perfect prediction detection, changing retained regressors
forvalues i=2/0`lc_n' { // reassemble hierarchical command line
local cmdline2 `cmdline2' || `lc_`i''
}
* re-estimate!
cap `=cond(`quietly', "", "noisily")' `cmdline' `from' `init' `iterate' `=cond("`cmd'"=="slogit","nocorner","")' constraints(`_constraints') `cmdline2'
local rc = _rc
constraint drop `_constraints'
if e(converged)==0 {
di as err "Could not impose null."
exit 430
}
if !`rc' {
mat `b0' = e(b)
_mkvec `b', update from(`b0', skip) first // mlexp doesn't handle from(..., skip)
cap `cmdline' `=cond(inlist("`cmd'","cmp","ml"),"init(`b')","from(`b')")' iterate(0) `cmdline2'
if _rc cap `cmdline' from(`b', skip) iterate(0) `cmdline2' // mprobit still needs skip even after _mkvec, from(..., skip)
local rc = _rc
mat drop `b0'
}
if `rc' {
di as err "Error imposing null. Perhaps {cmd:`cmd'} does not accept the {cmd:constraints()}, {cmd:from()}, and {cmd:iterate()} options, as needed."
exit `rc'
}
if !`quietly' {
mata st_numscalar("`t'", any(!diagonal(st_matrix("e(V)")) :& st_matrix("e(b)")'))
if `t' {
di _n as res _n "{res}Warning: Negative Hessian under null hypothesis not positive definite. Results may be unreliable." _n
di "This may indicate that the constrained optimum is far from the unconstrained one, i.e., that the hypothesis is incompatible with the data." _n
}
}
}
mat `V' = e(V`=cond("`e(V_modelbased)'"!="", "_modelbased", "")')
tempvar sc
cap predict double `sc'* if e(sample), score
if _rc {
if `reps'==0 & "`e(gradient)'"!="" { // hack: for score and Wald tests, generate observation-level "scores" with last entries=total gradient and rest 0, since only their sums matter
tempname gradient
mat `gradient' = e(gradient)
mata st_view(_boottestt=.,.,"`hold'"); st_local("n1", strofreal(max(selectindex(_boottestt))))
local scnames
qui forvalues i=1/`k' {
tempvar t
gen double `t' = 0
replace `t' = `gradient'[1,`i'] in `n1'
local scnames `scnames' `t'
}
}
else {
di as err "Could not generate scores from `e(cmd)' with {cmd:predict}."
exit _rc
}
}
else {
unab scnames: `sc'*
if `:word count `scnames'' == e(k_eq) { // generate score for each parameter from those for each equation
local colnames: colnames `b'
tokenize `:coleq `b''
local lasteq
local eq 0
local _sc
qui forvalues i=1/`k' {
if "``i''" != "`lasteq'" | "``i''"=="/" {
local ++eq
local lasteq: copy local `i'
}
local var: word `i' of `colnames'
if "`var'"=="_cons" | strmatch("`var'","*._cons") | "``i''"=="/" local _sc `_sc' `:word `eq' of `scnames'' // constant term or free parameter
else {
tempvar t
gen double `t' = `:word `eq' of `scnames'' * `var' if e(sample)
local _sc `_sc' `t'
}
}
local scnames: copy local _sc
}
}
if !`null' {
cap _estimates drop `hold'
_estimates hold `hold', restore // rename estimation as `hold' instead of generating new sample marker
}
}
else { // not ML
if `partial' mat `R' = `R', J(`df', `:word count `e(partial1)'' + e(partialcons), 0) // add blanks for partialled-out
mata st_matrix("`R'" , st_matrix("`R'" )[,_boottestp]) // put cols in standardized order & drop those for 0-variance vars
cap _estimates drop `hold'
_est hold `hold', restore
if `ar' {
local kEnd: word count `Xnames_endog'
local kEx : word count `Xnames_exog'
mata st_numscalar("`p'", rows(st_matrix("`R'"))!=`kEnd' | (`kEx'? any(st_matrix("`R'")[|.,.\.,`kEx'|]) : 0))
if `p' {
di as err "For Anderson-Rubin test, null hypothesis must constrain all the instrumented variables and no others."
exit 198
}
}
}
cap confirm var `hold'
if _rc {
di as err "Sample marker for the regression not found. Perhaps the data set was cleared and reconstructed."
exit _rc
}
if `margins' {
cap drop `touse'
gen byte `touse' = `hold' & `marginstouse'
local sample: copy local touse
}
else local sample: copy local hold
return local seed = cond("`seed'"!="", "`seed'", "`c(seed)'")
if !`julia' {
mata boottest_stata("`teststat'", "`df'", "`df_r'", "`p'", "`padj'", "`cimat'", "`plotmat'", "`peakmat'", `level', `ptolerance', ///
`ML', `LIML', 0`fuller', `K', `ar', `null', `scoreBS', `jk', "`weighttype'", "`ptype'", "`statistic'", ///
"`madjust'", `N_h0s', "`Xnames_exog'", "`Xnames_endog'", ///
"`Ynames'", `=cond(`ML', `" "`b'", "`V'" "', `" "","" "')', "`ZExclnames'", "`hold'", "`scnames'", `hasrobust', "`allclustvars'", `NBootClustVar', `NErrClustVar', ///
"`FEname'", `NFE', `FEdfadj', "`wtname'", "`wtype'", "`R1'", "`r1'", "`R'", "`r'", `reps', "`repsname'", "`repsFeasname'", `small', "`svmat'", "`dist'", ///
"`gridmin'", "`gridmax'", "`gridpoints'", `matsizegb', `quietly', "`b0'", "`V0'", "`svv'", "`NBootClustname'")
}
else {
foreach X in R R1 V {
if "``X''" != "" jl PutMatToMat ``X'', dest(`X')
else _jl: `X' = Matrix{Float`precision'}(undef,0,0);
}
foreach X in gridminvec gridmaxvec gridpointsvec {
if "``X''" != "" jl PutMatToMat ``X'', dest(`X')
else _jl: `X' = Matrix{Float`precision'}(undef,`=`df'',0);
}
foreach X in r r1 b {
if "``X''" != "" {
jl PutMatToMat ``X'', dest(`X')
_jl: `X' = vec(`X');
}
else _jl: `X' = Vector{Float`precision'}(undef,0);
}
foreach X in Xnames_exog Xnames_endog ZExclnames scnames allclustvars {
if "``X''" != "" jl PutVarsToMat ``X'' if `hold', dest(`X')
else _jl: `X' = Matrix{Float`precision'}(undef,0,0);
}
foreach X in Ynames wtname FEname {
if "``X''" != "" {
jl PutVarsToMat ``X'' if `hold', dest(`X')
_jl: `X' = dropdims(`X'; dims=2);
}
else _jl: `X' = Vector{Float`precision'}(undef,0);
}
_jl: FEname = iszero(`NFE') ? Vector{Int64}(undef,0) : Vector{Int64}(FEname);
_jl: allclustvars = iszero(`hasclust') ? Matrix{Int64}(undef,0,0) : Matrix{Int64}(allclustvars);
// jl: using JLD
// jl: @save "c:/users/drood/Downloads/tmp.jld" Ynames Xnames_exog Xnames_endog ZExclnames wtname allclustvars FEname scnames R r R1 r1 gridminvec gridmaxvec gridpointsvec b V
_jl: using Random; Random.seed!(rng, `=runiformint(0, 9007199254740992)'); // chain Stata rng to Julia rng
_jl: test = wildboottest!(Float`precision', R, r; resp=Ynames, predexog=Xnames_exog, predendog=Xnames_endog, inst=ZExclnames, ///
obswt=wtname, clustid=allclustvars, feid=FEname, scores=scnames, ///
R1, r1, ///
nbootclustvar=`NBootClustVar', nerrclustvar=`NErrClustVar', ///
issorted=true, ///
hetrobust=Bool(`hasrobust'), ///
fedfadj=`FEdfadj', ///
fweights = "`wtype'"=="fweight", ///
maxmatsize=`=0`matsizegb'', ///
ptype="`ptype'", ///
bootstrapc = "`statistic'"=="c", ///
liml=Bool(`LIML'), fuller=`=0`fuller'', `=cond(`K'<.,"kappa=`K',","")' ///
arubin=Bool(`ar'), small=Bool(`small'), ///
scorebs=Bool(`scoreBS'), ///
jk=Bool(`jk'), ///
reps=`reps', imposenull=Bool(`null'), level=`level'/100, ///
auxwttype=:`weighttype', ///
rtol=`ptolerance', ///
madjtype="`madjust'"=="" ? :none : :`madjust', nH0=`N_h0s', ///
ml=Bool(`ML'), beta=b, A=V, ///
gridmin=gridminvec, gridmax=gridmaxvec, gridpoints=gridpointsvec, ///
getdist = Bool("`svmat'"!=""), ///
getci = `level'<100 && "`cimat'" != "", getplot = "`plotmat'"!="", ///
getauxweights = "`svv'"!="", ///
rng=rng);
qui _jl: rand(rng, Int32) // chain Julia rng back to Stata to advance it replicably
set seed `r(ans)'
if "`plotmat'"!="" {
if `df'==1 {
jl GetMatFromMat `plotmat', source([test.plot[:X][1] test.plot[:p]])
jl GetMatFromMat `peakmat', source([test.peak[:X][1] test.peak[:p]])
}
else jl GetMatFromMat `plotmat', source([vcat(vec([[x y] for x in test.plot[:X][1], y in test.plot[:X][2]])...) test.plot[:p]])
}
if `level'<100 & "`cimat'" != "" jl GetMatFromMat `cimat', source(test.ci)
_jl: SF_scal_save("`teststat'", test.stat);
_jl: SF_scal_save("`df'", test.dof);
_jl: SF_scal_save("`df_r'", test.dof_r);
_jl: SF_scal_save("`p'", test.p);
_jl: SF_scal_save("`padj'", test.padj);
_jl: SF_scal_save("`repsname'", test.reps);
_jl: SF_scal_save("`repsFeasname'", test.repsfeas);
_jl: SF_scal_save("`NBootClustname'", test.nbootclust);
jl GetMatFromMat `b0', source(test.b)
jl GetMatFromMat `V0', source(test.V)
if "`dist'"!="" jl GetMatFromMat `dist', source(test.`=cond("`svmat'"=="t", "dist", "numerdist")')
if "`svv'" !="" jl GetMatFromMat `svv', source(test.auxweights)
}
_estimates unhold `hold'
if !`quietly' & 2^`NBootClustname' < `reps' & inlist("`weighttype'", "rademacher", "mammen") {
di _n "Warning: with " `NBootClustname' " bootstrap clusters, the number of replications, `reps', exceeds the universe of " strproper("`weighttype'") " draws, 2^"`NBootClustname' " = " 2^`NBootClustname' ". " _c
if "`weighttype'"=="rademacher" di "Sampling each once." _c
di _n "Consider Webb weights instead, using {cmd:weight(webb)}."
}
local reps = `repsname' // in case reduced to 2^G
if !`quietly' & `reps' & (`NBootClustname'>12 | "`weighttype'" != "rademacher") & floor(`level'/100 * (`reps'+1)) != `level'/100 * (`reps'+1) {
di _n "Note: The bootstrap usually performs best when the confidence level (here, `level'%)" _n " times the number of replications plus 1 (`reps'+1=" `reps'+1 ") is an integer."
}
if !`ML' & `reps' & `NErrClustVar'>1 & `teststat'<. {
return scalar repsFeas = `repsFeasname'
if `repsFeasname' < `reps' di _n "Warning: " `reps' - `repsFeasname' " replications returned an infeasible test statistic and were deleted from the bootstrap distribution."
}
if `N_h0s'>1 local _h _`h'
if `df'==1 {
local tz = cond(`small', "t", "z")
return scalar `tz'`_h' = `teststat'
}
di
if `reps' {
di as txt strproper("`boottype'") " bootstrap-`statistic', null " cond(0`null', "", "not ") "imposed, " _c
if `jk' di "jackknifed residuals, " _c
di as txt `reps' as txt " replications, " _c
}
di as txt cond(`ar', "Anderson-Rubin ", "") cond(!`reps' & `null' & "`boottype'"=="score", "Rao score (Lagrange multiplier)", "Wald") " test" _c
if "`cluster'"!="" di ", clustering by " as res "`cluster'" _c
if ("`bootcluster'"!="" | `NErrClustVar' > 1) & `reps' di as txt ", bootstrap clustering by " as res "`bootcluster'" _c
if `reps' di as txt ", " strproper("`weighttype'") " weights" _c
di as txt ":"
if `margins' di `" `:word `h' of "`marginsnames'""'
else {
foreach c in `h0_`h'' {
di " " _c
if !0`anythingh0' di as txt %4.0g `c' ": " _c
di as res "`:constraint `c''"
}
}
if `df'==1 {
local line1 = cond(`small', "t(`=`df_r'')", "z")
di _n as txt _col(`=33-strlen("`line1'")') "`line1' = " as res %10.4f `teststat' _n _col(`=33-strlen("``ptype''")') as txt "``ptype'' = " as res %10.4f `p'
local `teststat' = `teststat' * `teststat'
}
else {
if `small' {
local line1 F(`=`df'', `=`df_r'')
local line2 Prob > F
}
else {
local line1 chi2(`=`df'')
local line2 Prob > chi2
}
di _n as txt _col(`=27-strlen("`line1'")') "`line1' = " as res %10.4f `teststat' _n as txt _col(`=27-strlen("`line2'")') "`line2' = " as res %10.4f `p'
}
if "`madjust'" != "" {
di _col(`=strlen("bonferroni")-strlen("`madjust'")+3') as txt strproper("`madjust'") "-adjusted prob = " as res %10.4f `padj'
return scalar padj`_h' = `padj'
}
if `NErrClustVar' > 1 {
cap mat `t' = syminv(`V0')
cap noi if diag0cnt(`t') di _n "Test statistic could not be computed. The multiway-clustered variance estimate " cond(`df'==1, "", "matrix ") "is not positive" cond(`df'==1, "", "-definite") "."
}
cap mat colnames `b0' = `h0text'
cap mat colnames `V0' = `h0text'
cap mat rownames `V0' = `h0text'
cap return matrix b`_h' = `b0'
cap return matrix V`_h' = `V0'
cap confirm mat `plotmat'
if _rc == 0 {
tempvar X1 Y _plotmat
mat `_plotmat' = `plotmat'
return matrix plot`_h' = `plotmat'
mat `plotmat' = `_plotmat'
if "`graph'"=="" {
cap confirm matrix `peakmat'
if !_rc {
mata st_matrix("`plotmat'", st_matrix("`plotmat'") \ st_matrix("`peakmat'"))
return matrix peak`_h' = `peakmat'
}