-
Notifications
You must be signed in to change notification settings - Fork 3
/
lassoutils.ado
6154 lines (5617 loc) · 197 KB
/
lassoutils.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
*! lassoutils 1.2.06 5jan2024
*! lassopack package 1.4.3
*! authors aa/cbh/ms
* Adapted/expanded from lassoShooting version 12, cbh 25/09/2012
* mm_quantile from moremata package by Ben Jann:
* version 1.0.8 20dec2007 Ben Jann
* Notes:
* Partialling out, temp vars, factor vars, FE transform all handled
* by calling programs.
* names_o and names_t (via options) are original and (possibly) temp names of Xs
* lassoutils branches to internal _rlasso, _lassopath, _fe, _unpartial, _partial, _std
* current cluster code is memory-intensive since it works with a full n x p matrix instead of cluster-by-cluster
* Updates (release date):
* 1.0.05 (30jan2018)
* First public release.
* Added seed(.) option to rlasso to control rnd # seed for xdep & sup-score.
* Fixed up return code for lassoutils (method, alpha).
* Promoted to required version 13 or higher.
* Introduced centerpartial(.) Mata program for use with rlasso; returns X centered and with Xnp partialled-out.
* Separate fields in datastruct: sdvec, sdvecpnp, sdvecp, sdvecnp. Latter two conformable with Xp and Xnp.
* Added dots option for simulations (supscore, xdep).
* Recoding relating to different treatment of cross-validation.
* Changes to _std, _fe, _partial relating to holdout vs full sample.
* Recoding of cons flag; now also dmflag to indcate zero-mean data.
* 1.0.06 (10feb2018)
* Misc Mata speed tweaks following advice at http://scorreia.com/blog/2016/10/06/mata-tips.html,
* e.g., evaluating for limits before loop; referring to vector elements with a single subscript.
* Rewrite of FE transform to use Sergio Correia's FTOOLS package (if installed).
* 1.0.07 (17feb2018)
* Bug fix related to ftools - leaves matastrict set to on after compilation, causing rest of lassoutils
* to fail to load. Fix is to reset matastrict to what it was before calling ftools,compile.
* 1.0.08 (5apr2018)
* Changed the default maximum lambda for the elastic net: it was 2*max(abs(X'y)),
* and is now 2*max(abs(X'y))/max(0.001,alpha); see Friedman et al (2010, J of Stats Software).
* Added Mata programs getInfoCriteria() and getMinIC(). getInfoCriteria() calculates information criteria
* along with RSS and R-squared. getMinIC() obtains minimum IC values and corresponding lambdas.
* Degrees of freedom calculation was added to DoLassoPath() and DoSqrtLassoPath().
* Misc changes to the initial beta (= Ridge estimate), which in some cases didn't account for
* penalty loadings.
* Added lglmnet option to facilitate comparison with glmnet (was dysfunctional 'ladjust' option).
* Restructured calc of lambda with cluster to match JBES paper; prev used sqrt(nclust) rather than sqrt(N),
* now always 2c*sqrt(N)*invnormal(1-(gamma/log(N))/(2*p))). Fixed bug in calc of lambda for sqrt-lasso with cluster.
* Definition of xdep lambda + cluster also changed slightly (by factor of sqrt(nclust/(nclust-1)).
* Undocumented nclust1 option mimics behavior of CBH lassoCluster.ado; uses (nclust-1)*T rather than nclust*T.
* 1.0.14 (04sep2018)
* Bug with initial (ridge) beta if XX+lambda2/2*diag(Ups2) singular, caused by using lusolve.
* Fixed with subroutine luqrsolve(.). Uses lusolve by default; if singular, then uses qrsolve.
* Dropped unpenalized variables now has more informative error message.
* Modified code to accommodate cluster variables that are strings.
* Added weighting subroutine (data assumed to be preweighted).
* Consistent use of d.dmflag and d.cons:
* dmflag=1 => treat data as zero mean
* cons=1 => constant in model to be recovered by estimation code; also used in variable counts
* dmflag=1 => cons=0 but cons=0 /=> dmflag=1 because model may not have a constant
* Fixed wrong formula for R-squared in getInfoCriteria(). RSQ =1:-RSS:/TSS (was ESS:/TSS).
* Fixed initial resids routine for rlasso for model with no constant.
* Fixed bug in sup score test + prestd. Fixed bug that didn't allow sup score test + pnotpen(.).
* Added separate storage of prestandardization SDs on data struct as prestdx and prestdy.
* Added "ebicgamma" option. Default choice is now ebicgamma=1-1/(2*kappa) and p=n^kappa.
* See Chen & Chen (2008, Sec. 5, p. 768).
* 1.1.01 (08nov2018)
* Rewrite of supscore code. Now reports sqrt(n)*[] rather than n*[] version as per CCK.
* Bug fix in cluster version (was using nclust*[] rather than n*[]).
* Now handles weighted/demeaned data, center option correctly.
* Rewrite of lambdacalc code. Shares subroutines with supscore code.
* Rewrite of xdep and lasso/sqrt-lasso weights code. Shares subroutines with supscore code.
* sqrt-lasso xdep simulation now normalizes by empirical SD of simulated N(0,1).
* Bug fix in xdep gamma; distribution was based on max(Xe) instead of max(abs(Xe));
* effect was appx equivalent to using a gamma 1/2 times the value specified; now fixed.
* Rewrite of cluster weights code to use panelsum(.) and avoid large temp matrices.
* Added saved value of minimized objective function (standardized).
* Initial weights for sqrt-lasso + inid data as per BCW 2014 p. 769 = colmax(abs(X))
* now available as undocumented maxabsx option. For cluster, use colmax(panelmean(abs(X)).
* maxabsx => default maxupsiter=3 in line with BCW 2014 recommendation to iterate.
* Bug fix in reporting of sqrt-lasso with vverbose (wasn't reporting correct value of obj fn).
* Bug fix in sqrt-lasso with nocons (was demeaning despite no constant present).
* Bug fix in rlasso verbose reporting - wasn't reporting selected variables.
* Removed factor of 0.5 from initial penalty (loadings) for lasso + post-lasso OLS resids.
* Added undocumented c0(.) option to replicate previous factor of 0.5 for initial penalty loadings.
* Added undocumented sigma(.) option for std lasso to override estimation of initial sigma.
* Added undocumented ssiid option for supscore test to override use of mult bootstrap in iid case.
* 1.1.02 (2jan2018)
* Saved value of minimized objective function is now in y units (not standardized).
* Adaptive lasso: added check for zeros in e(b). Use univariate OLS instead.
* 1.1.03 (13jan2019)
* Updated options and terminology from Ups to Psi (penalty loadings).
* gamma option now controls overall level, = 0.1/log(N) or 0.1/log(N_clust) by default
* gammad option now undocumented
* 1.1.04 (14oct2019)
* Bug fix - lambda0(.) option was being ignored for standard lasso.
* Algorithm update for rlasso - now maxupsiter=1 will exit with the initial penalty loadings,
* initial lambda, and lasso estimates based on this.
* NB: Also specifying lambda0(.) will control the initial lambda.
* 1.2.01 (29july2020)
* Substantial rewrite and simplification of code, misc bug fixes, new features.
* HAC/AC rlasso added. Imports code from livreg2 with modifications.
* 2-way cluster-robust added. rlasso now also returns r2.
* DoSqrtLasso(.) now incorporated into DoLasso(.); similarly for DoSqrtLassoPath(.).
* Shooting algo convergence now in standardized metric
* DoLassoPath(.) lambda iteration report; small coefs not set to zero if rige; compares to std coeffs;
* saves untruncated lambda list; dof and model dimension fixes (counts all partialled out or none);
* exits path early if dev criteria met; override and obtain old behavior with nodevcrit option.
* Bug fix in lassopath with stdcof+supplied lamba.
* lassopath now exits if lminratio is illegal.
* lassopath saves shat0 (shat+partialled out) and lambdamat0 (untruncated lambdamat).
* Bug fixes in info criteria code (check for constant, handling of sqrt-lasso, demeaning of y).
* Effective degrees of freedom bug fix (accounts for constant and partialled-out).
* Adalasso fixes: initialization of zerosinada; scaling of adaptive weights; varnames checks;
* initialization of Psi (default 1e-9 rather than exact 0); adaptive now matches adaptive+prestd;
* abs(.) precedes raising to theta power; verbose behavior
* Unpartial bug fix (could crash if no selected regressors).
* Solver option allowed for partialling out; recoded to be string rather than numeric.
* Used for both partialling out by calling program and for rlasso partialling.
* Default solver for partialling is now qrxx = QR+quadcross (was SVD+QR).
* Warning issued if collinearities encountered when partialling out.
* Bug fix in case all ICs were missing.
* nodevcrit option (uses full lambda list irrespective of size of deviance in path)
* 1.2.02 (4aug2020)
* Bug fix in adaptive lasso relating to omitted vars (would get a missing penalty loading).
* Reverted to adaptive lasso behavior where any omitted vars trigger univariate OLS.
* Bug fix in lmax+prestd option (lmax should be rescaled).
* Reporting bug fix with vverbose option and elastic net - wrong obj fn. (Reporting only.)
* 1.2.03 (27sept2020)
* Added lglmnet option for _lassopath.
* Recoded algorithms to drop use of lambda2 and psi2.
* Removed * for options so illegals caught.
* Bug fix in value returned for ridge objective function (was just rmse, excluded penalty).
* Bug fix in starting ridge coefs for rlasso.
* Added check for negative penalty loadings.
* Bug fix for user-supplied penalty loadings and adaptive lasso with alpha<1.
* Calculation of intercepts moved to return results subroutines.
* Adaptive lasso recoded in Mata.
* Adaptive elastic net now correctly uses adaptive weights for L1 norm only.
* Standardized and unstandardized coefficients now returned for lassopath; stdcoef option unneeded.
* Standardized lambdas, norms, ICs, etc. also returned.
* lassopath code takes stdall option - indicates that provided lambdas are in the standardized metric
* Fixed bug in reporting of value of maximized obj fn for elastic net/ridge (missing 1/2 on L2 norm).
* EBIC now excludes omitted/base variables when calculating model size p.
* 1.2.04 Bug fix to SD calculation for special case of lglmnet with unit loadings (=not prestandardized)
* 1.2.05 Bug fix in special case of only unpenalized regressors (returned scalar psinegs was not initialized to missing)
* 1.2.06 (5jan2024)
* Misc updates to support use of Python/sklearn.
program lassoutils, rclass sortpreserve
version 13
syntax [anything] , /// anything is actually a varlist but this is more flexible - see below
[ ///
rlasso /// branches to _rlasso
path /// branches to _lassopath
unpartial /// branches to _unpartial
fe(string) /// branches to _fe
std /// branches to _std
wvar(string) /// branches to _wt
partialflag(int 0) /// branches to _partial if =1
partial(string) /// used by _partial and _unpartial
tvarlist(string) /// used by _partial, _fe, _std; optional list of temp vars to fill
tpvarlist(string) /// used by _partial; optional list of temp partial vars to file
///
ALPha(real 1) /// elastic net parameter
SQRT /// square-root-lasso
ols ///
adaptive ///
///
/// verbose options will be transformed to 0/1/2 before passing to subs
VERbose /// shows additional detail
VVERbose /// even more detail
* ///
]
*** `anything' is actually a varlist but this is more flexible.
// avoids problems with uninitialized vars, unwanted parsing of factor vars, etc.
// ... so immediately rename.
local varlist `anything'
*
** verbose
if ("`vverbose'"!="") {
local verbose=2 // much detail
}
else if ("`verbose'"!="") {
local verbose=1 // some additional detail
}
else {
local verbose=0 // no additional output
}
*
************* SUBROUTINES: shoot, path, CV **************************************
if "`rlasso'" ~= "" {
_rlasso `varlist', /// branch to _rlasso
verbose(`verbose') ///
`sqrt' ///
`options' // no penloads etc.
}
else if "`path'" != "" {
_lassopath `varlist', /// branch to _ lassopath
verbose(`verbose') ///
`sqrt' ///
alpha(`alpha') ///
`adaptive' ///
`ols' ///
`options'
}
else if "`unpartial'" ~= "" {
_unpartial `varlist', /// branch to _unpartial
partial(`partial') ///
wvar(`wvar') ///
`options'
}
else if `partialflag' {
_partial `varlist', /// branch to _partial
partial(`partial') ///
tvarlist(`tvarlist') ///
wvar(`wvar') ///
`options'
}
else if "`fe'" ~= "" {
_fe `varlist', /// branch to _fe
tvarlist(`tvarlist') ///
fe(`fe') ///
wvar(`wvar') ///
`options'
}
else if "`std'" ~= "" {
_std `varlist', /// branch to _std
tvarlist(`tvarlist') ///
`options'
}
else if "`wvar'" ~= "" {
_wt `varlist', /// branch to _wt
tvarlist(`tvarlist') ///
wvar(`wvar') ///
`options'
}
else {
di as err "internal lassoutils error"
exit 499
}
*** return
// so subroutine r(.) detail passed on
// can add additional items to r(.) here
return add
// if estimation, return method etc.
if "`rlasso'`path'`cv'" ~= "" {
return scalar alpha = `alpha'
return scalar sqrt = ("`sqrt'"!="")
if ("`sqrt'" ~= "") {
return local method "sqrt-lasso"
}
else if `alpha'==1 {
return local method "lasso"
}
else if `alpha'==0 {
return local method "ridge"
}
else if (`alpha'<1) & (`alpha'>0) {
return local method "elastic net"
}
else {
di as err "internal lassoutils error - alpha outside allowed range"
exit 499
}
}
end
// Subroutine for BCH _rlasso
program define _rlasso, rclass sortpreserve
version 13
syntax anything , /// anything is actually a varlist but this is more flexible - see below
touse(varlist numeric max=1) /// required `touse' variable (so marksample unneeded)
[ ///
verbose(int 0) ///
///
ROBust /// penalty level/loadings allow for heterosk. [homoskedasticity is the default]
CLuster(varlist max=2) /// penalty level/loadings allow for within-panel dependence & heterosk.
bw(int 0) /// HAC bandwidth
kernel(string) ///
tindex(varlist max=1) ///
bsize(int 5) /// default block size for the HAC/AC multiplier bootstrap
maq /// specific to HAC/AC with truncated kernel
XDEPendent /// penalty level is estimated depending on X (default: No)
numsim(integer 5000) /// number of simulations with xdep (default=5,000)
///
TOLOpt(real 1e-10) /// was originally 1e-5 but 1e-10 gives good equiv between xstd and stdloadings
TOLPsi(real 1e-4) ///
TOLZero(real 1e-4) ///
MAXIter(integer 10000) /// max number of iterations in estimating lasso
MAXPSIIter(int 2) /// max number of lasso-based iterations in est penalty loadings; default=2
CORRNumber(int 5) /// number of regressors used in InitialResiduals(); default=5
maxabsx /// initial loadings for sqrt-lasso as per BCW 2014 p. 769 = colmax(abs(X))
///
LASSOPSI /// use lasso residuals to estimate penalty loadings (post-lasso is default)
///
c(real 1.1) /// "c" in lambda function; default = 1.1
c0(real 0) /// when iterating, initial value of c in lambda function; default set below
gamma(real 0) /// "gamma" in lambda function (default 0.1/log(N) or 0.1/log(N_clust))
gammad(real 0) /// undocumented "gamma" denominator option (default log(N) or log(N_clust)
lambda0(real 0) /// optional user-supplied lambda0
LALTernative /// use alternative, less sharp lambda0 formula
PMINus(int 0) /// dim(X) adjustment in lambda0 formula
nclust0 /// no longer in use as of 1.0.08; replaced by nclust1
nclust1 /// use (nclust-1)*T instead of nclust*T in cluster-lasso
center /// center x_i*e_i or x_ij*e_ij
sigma(real 0) /// user-specified sigma for rlasso
///
supscore ///
ssgamma(real 0.05) /// default gamma for sup-score test
ssnumsim(integer 500) ///
testonly ///
ssiid /// undocumented option - override use of mult bootstrap in iid case
///
seed(real -1) /// set random # seed; relevant for xdep and supscore
dots /// display dots for xdep and supscore repetitions
///
xnames_o(string) /// original names in varXmodel if tempvars used (string so varlist isn't processed)
xnames_t(string) /// temp names in varXmodel
notpen_o(string) /// used in checking
notpen_t(string) ///
consflag(int 0) /// =0 if no cons or already partialled out
dofminus(int 0) /// lost degrees of freedom from FEs
sdofminus(int 0) /// lost degrees of freedom from partialling out ("s"=small; includes the constant)
dmflag(int 0) /// data have been demeaned or should be treated as demeaned
///
stdy(string) ///
stdx(string) ///
///
SQRT /// square-root lasso
///
/// LEGACY OPTIONS
TOLUps(real 1e-4) ///
MAXUPSIter(int 2) /// max number of lasso-based iterations in est penalty loadings; default=2
LASSOUPS /// use lasso residuals to estimate penalty loadings (post-lasso is default)
psolver(string) ///
]
*** `anything' is actually a varlist but this is more flexible.
// avoids problems with uninitialized vars, unwanted parsing of factor vars, etc.
// ... so immediately rename.
local varlist `anything'
*
*** count number of obs
_nobs `touse'
local nobs = r(N)
*
*** cluster
// create numeric clustvar in case cluster variable is string
// sort needed if cluster option provided; will need to re-sort afterwards
local clustercount : word count `cluster'
if `clustercount' {
local sortedby : sortedby
}
if `clustercount'==1 {
tempvar clustvar
qui egen `clustvar' = group(`cluster') if `touse'
sort `clustvar'
local allclust `clustvar'
}
else if `clustercount'==2 {
tokenize `cluster'
tempvar clustvar clustvar2 clustvar3
qui egen `clustvar' = group(`1') if `touse'
qui egen `clustvar2' = group(`2') if `touse'
qui egen `clustvar3' = group(`1' `2') if `touse'
sort `clustvar3' `clustvar'
local allclust `clustvar' `clustvar2' `clustvar3'
}
*
*** HAC/AC
if `bw'~=0 {
capture tsset
if "`r(timevar)'" == "" {
di as err "must tsset data and specify timevar"
exit 5
}
local tdelta = `r(tdelta)'
// check if valid bandwidth
sum `r(timevar)' if `touse', meanonly
local T = r(max)-r(min) + 1
if `bw' > (`T'-1)/`tdelta' {
di as err "invalid bandwidth in option bw() - cannot exceed timespan of data"
exit 198
}
// check if valid kernel, standardize kernel name, also set spectral flag
mata: s_vkernel("`kernel'", "`bw'")
local kernel `r(kernel)'
local spectral = `r(spectral)'
local maqflag = ("`maq'"~="" & "`kernel'"=="Truncated")
}
else {
// bw=0 means ignore kernel option and set spectral to default of 0
local tdelta = 0
local kernel
local spectral = 0
local maqflag = 0
}
*
// incompatible options
if `clustercount' & `bw' {
di as err "incompatible options - cluster(.) and bw(.)"
exit 198
}
*** define various parameters/locals/varlists
local varY `varlist' // now is only element in varlist
local varXmodel `xnames_t' // including notpen if any; may include tempnames
local pen_t : list xnames_t - notpen_t // list of penalized regressors
local p0 : word count `varXmodel' // #regressors incl notpen but excl constant
local p = `p0' - `pminus' // p defined by BCH
local p0 = `p0' + `consflag' // #regressors but now incl constant
*
*** define various flags
local clustflag =("`cluster'"!="") // =1 if cluster
local hetero =("`robust'"!="") & ~`clustflag' // =1 if het-robust but NOT cluster
local sqrtflag =("`sqrt'"!="")
local xdep =("`xdependent'"!="")
local lassopsiflag =("`lassopsi'`lassoups'"!="") // lassoups is the equivalent legacy option
local lalternative =("`lalternative'"!="") // =1 if alternative, less sharp lambda0 should be used
local nclust1flag =("`nclust1'"!="")
local center =("`center'"!="")
local supscoreflag =("`supscore'`testonly'"!="")
local testonlyflag =("`testonly'"!="")
local dotsflag =("`dots'"!="")
local prestdflag =("`stdy'"!="")
local ssiidflag =("`ssiid'"!="")
local maxabsxflag =("`maxabsx'"!="")
*
*** defaults
if `maxabsxflag' {
// implement BCW 2014 p. 769 recommendation for sqrt lasso initial weights = colmax(abs(X))
local corrnumber -1
local maxpsiiter 3 // since BCW say to iterate
}
if `c0'==0 {
// allow initial value of c in iterated penalties to vary; matches treatment in CBH and hdm code
local c0 `c'
}
*
*** legacy options
if `tolups' ~= 1e-4 {
local tolpsi `tolups'
}
if `maxupsiter' ~= 2 {
local maxpsiiter `maxupsiter'
}
*
tempname lambda slambda rmse rmseOLS objfn r2 // note lambda0 already defined as a macro
tempname b bOLS sb sbOLS Psi sPsi ePsi stdvec
tempname bAll bAllOLS
tempname supscoremat CCK_ss CCK_p CCK_cv CCK_gamma
tempname psinegs
// initialize so that returns don't crash in case values not set
local k =0 // used as flag to indicate betas are non-missing
local N =.
local N_clust =.
local s =.
local s0 =.
local niter =.
local npsiiter =.
scalar `lambda' =.
scalar `slambda' =.
scalar `rmse' =.
scalar `rmseOLS' =.
scalar `r2' =.
scalar `objfn' =.
mat `b' =.
mat `bAll' =.
mat `bOLS' =.
mat `bAllOLS' =.
mat `sb' =.
mat `sbOLS' =.
mat `Psi' =.
mat `sPsi' =.
mat `ePsi' =.
mat `stdvec' =.
mat `supscoremat' =.
scalar `CCK_ss' =.
scalar `CCK_p' =.
scalar `CCK_cv' =.
scalar `CCK_gamma' =.
scalar `psinegs' =.
if `p' & ~`testonlyflag' { // there are penalized model variables; estimate
*** BCH rlasso
mata: EstimateRLasso( ///
"`varY'", ///
"`varXmodel'", ///
"`xnames_o'", ///
"`pen_t'", ///
"`notpen_t'", /// #5
"`touse'", ///
"`allclust'", ///
"`stdy'", ///
"`stdx'", ///
`sqrtflag', /// #10
`hetero', ///
`bw', ///
"`kernel'", ///
`maqflag', ///
`spectral', /// #15
"`tindex'", ///
`tdelta', ///
`bsize', ///
`xdep', ///
`numsim', /// #20
`lassopsiflag', ///
`tolopt', ///
`maxiter', ///
`tolzero', ///
`maxpsiiter', /// #25
`tolpsi', ///
`verbose', ///
`c', ///
`c0', ///
`gamma', /// #30
`gammad', ///
`lambda0', ///
`lalternative', ///
`corrnumber', ///
`pminus', /// #35
`nclust1flag', ///
`center', ///
`sigma', ///
`supscoreflag', ///
`ssnumsim', /// #40
`ssgamma', ///
`seed', ///
`dotsflag', ///
`consflag', ///
`dmflag', /// #45
`dofminus', ///
`sdofminus', ///
`prestdflag', ///
"`psolver'")
*** Via ReturnResults(.)
// coefs are returned as column vectors
// convert to row vectors (Stata standard)
mat `b' =r(b)' // in original units
mat `bOLS' =r(bOLS)'
mat `sb' =r(sb)' // in standardized units
mat `sbOLS' =r(sbOLS)'
mat `bAll' =r(bAll)'
mat `bAllOLS' =r(bAllOLS)'
mat `Psi' =r(Psi)
mat `sPsi' =r(sPsi)
mat `ePsi' =r(ePsi)
mat `stdvec' =r(stdvecp) // stdvec for penalized vars only
local selected0 `r(sel)' // selected variables INCL NOTPEN, EXCL CONSTANT
local s0 =r(s) // number of selected vars INCL NOTPEN, EXCL CONSTANT; may be =0
local k =r(k) // number of all vars in estimated parameter vector INC CONSTANT; may be =0
local niter =r(niter)
local npsiiter =r(npsiiter)
local N =r(N)
local N_clust =r(N_clust)
local N_clust1 =r(N_clust1)
local N_clust2 =r(N_clust2)
scalar `lambda' =r(lambda) // relates to depvar in original units
scalar `slambda' =r(slambda) // relates to standardized depvar
scalar `rmse' =r(rmse) // lasso rmse
scalar `rmseOLS' =r(rmsePL) // post-lasso rmse
scalar `r2' =r(r2)
scalar `objfn' =r(objfn) // minimized objective function
if `supscoreflag' {
mat `supscoremat' =r(supscore) // sup-score row vector of results
scalar `CCK_ss' =`supscoremat'[1,colnumb(`supscoremat',"CCK_ss")]
scalar `CCK_p' =`supscoremat'[1,colnumb(`supscoremat',"CCK_p")]
scalar `CCK_cv' =`supscoremat'[1,colnumb(`supscoremat',"CCK_cv")]
scalar `CCK_gamma' =`supscoremat'[1,colnumb(`supscoremat',"CCK_gamma")]
}
// these may overwrite existing locals
tempname lambda0 c gamma gammad
scalar `lambda0' =r(lambda0) // BCH definition of lambda; overwrites existing/default macro
scalar `c' =r(c)
scalar `gamma' =r(gamma)
scalar `gammad' =r(gammad)
// for HAC or 2-way cluster lasso (neg loadings possible)
scalar psinegs =r(psinegs)
local psinegvars `r(psinegvars)'
*
}
else if ~`testonlyflag' { // there are no penalized model vars; just do OLS
if `p0' { // there are unpenalized vars and/or constant
di as text "warning: no penalized variables; estimates are OLS"
}
if ~`consflag' {
local noconstant noconstant
}
else {
local consname _cons
}
qui regress `varY' `varXmodel' if `touse', `noconstant'
mat `b' =e(b)
local k =colsof(`b')
foreach m in `b' `bOLS' `bAll' `bAllOLS' {
mat `m' =e(b)
mat colnames `m' =`xnames_o' `consname'
}
scalar `rmse' =e(rmse)
scalar `rmseOLS' =e(rmse)
scalar `r2' =e(r2)
local selected0 `xnames_o'
local N =e(N)
local N_clust =.
local N_clust1 =.
local N_clust2 =.
}
else if `testonlyflag' { // supscore test only
mata: EstimateSupScore( ///
"`varY'", ///
"`varXmodel'", ///
"`xnames_o'", ///
"`pen_t'", ///
"`notpen_t'", /// #5
"`touse'", ///
"`allclust'", ///
"`stdy'", ///
"`stdx'", ///
`sqrtflag', /// #10
`hetero', ///
`bw', ///
"`kernel'", ///
`maqflag', ///
`spectral', /// #15
"`tindex'", ///
`tdelta', ///
`bsize', ///
`verbose', ///
`ssnumsim', /// #20
`ssiidflag', ///
`c', ///
`nclust1flag', ///
`center', ///
`ssgamma', /// #25
`pminus', ///
`seed', ///
`dotsflag', ///
`consflag', ///
`dmflag', /// #30
`dofminus', ///
`sdofminus', ///
`prestdflag', ///
"`psolver'")
mat `supscoremat' =r(supscore)
scalar `CCK_ss' =`supscoremat'[1,colnumb(`supscoremat',"CCK_ss")]
scalar `CCK_p' =`supscoremat'[1,colnumb(`supscoremat',"CCK_p")]
scalar `CCK_cv' =`supscoremat'[1,colnumb(`supscoremat',"CCK_cv")]
scalar `CCK_gamma' =`supscoremat'[1,colnumb(`supscoremat',"CCK_gamma")]
local N =r(N)
local N_clust =r(N_clust)
local N_clust1 =r(N_clust1)
local N_clust2 =r(N_clust2)
// for HAC or 2-way cluster lasso (neg loadings possible)
scalar psinegs =r(psinegs)
local psinegvars `r(psinegvars)'
}
else { // shouldn't reach here
di as err "internal _rlasso error"
exit 499
}
**************** Estimation complete ********************
// cluster required sorting, so need to restore original sort if was originally sorted
if `clustercount' & "`sortedby'"~="" {
sort `sortedby'
}
*** check notpen is in selected0
if ~`testonlyflag' {
fvstrip `selected0' // fvstrip removes b/n/o prefixes (shouldn't need dropomit)
local list1 `r(varlist)'
fvstrip `notpen_o', dropomit // use dropomit option here
local list2 `r(varlist)'
local missing : list list2 - list1
local missing_ct : word count `missing'
if `missing_ct' {
di as err "internal _rlasso error - unpenalized `missing' missing from selected vars"
di as err "may be caused by collinearities in unpenalized variables or by tolerances"
di as err "set tolzero(.) or other tolerances smaller or use partial(.) option"
exit 499
}
}
*
*** conventions
// k = number of selected/notpen INCLUDING constant; k>0 means estimated coef vector is non-empty
// s0 = number of selected/notpen EXCLUDING constant
// s = number of selected (penalized)
// p0 = number of all variables in model INCLUDING constant
// selected0 = varlist of selected and unpenalized EXCLUDING constant; has s0 elements
// selected = varlist of selected; has s elements
// cons = 1 if constant, 0 if no constant
// notpen_ct = number of unpenalized variables in the model EXCLUDING CONSTANT
// coef vectors b, bOLS, sb, sbOLS have k elements
// coef vectors bAll, bAllOLS have p0 elements
// Psi, sPsi have p0-cons elements
// stdvec has p0 - notpen_ct - cons elements (=#penalized)
// Note that full set of regressors can including omitteds, factor base categories, etc.
*
*** fix colnames of beta vectors to include omitted "o." notation
// includes post-lasso vector in case of collinearity => a selected variable was omitted
// trick is to use _ms_findomitted utility but give it
// diag(bAll) as vcv matrix where it looks for zeros
// also build in fv info
tempname tempvmat
if `k' & ~`testonlyflag' { // if any vars in est coef vector (k can be zero)
mat `tempvmat' = diag(`bOLS')
_ms_findomitted `bOLS' `tempvmat'
_ms_build_info `bOLS' if `touse'
}
if `p0' & ~`testonlyflag' { // if any variables in model (p0 can be zero)
mat `tempvmat' = diag(`bAll')
_ms_findomitted `bAll' `tempvmat'
_ms_build_info `bAll' if `touse'
mat `tempvmat' = diag(`bAllOLS')
_ms_findomitted `bAllOLS' `tempvmat'
_ms_build_info `bAllOLS' if `touse'
}
*
*** manipulate variable counts and lists
// selected0 and s0 are selected+notpen (excluding constant)
// selected and s are selected only
local notpen_ct : word count `notpen_o' // number of notpen EXCL CONSTANT
local selected : list selected0 - notpen_o // selected now has only selected/penalized variables
local s : word count `selected' // number of selected/penalized vars EXCL NOTPEN/CONSTANT
*
*** error checks
// check col vector of b (selected) vs. k
local col_ct =colsof(`b')
if `col_ct'==1 & el(`b',1,1)==. & ~`testonlyflag' {
// coef vector is empty so k should be zero
if `k'~=0 {
di as err "internal _rlasso error - r(k)=`k' does not match number of selected vars/coefs=`col_ct'"
exit 499
}
}
else if `k'~=`col_ct' & ~`testonlyflag' {
// coef vector not empty so k should match col count
di as err "internal _rlasso error - r(k)=`k' does not match number of selected vars/coefs=`col_ct'"
exit 499
}
// check col vector of bAll vs. p0
local col_ct =colsof(`bAll')
if `p0'~=`col_ct' & ~`testonlyflag' {
// full coef vector not empty so p0 should match col count
di as err "internal _rlasso error - p0=`p0' does not match number of model vars/coefs=`col_ct'"
exit 499
}
*
*** return results
// coef vectors are row vectors (Stata standard)
return matrix beta =`b'
return matrix betaOLS =`bOLS'
return matrix sbeta =`sb'
return matrix sbetaOLS =`sbOLS'
return matrix betaAll =`bAll'
return matrix betaAllOLS=`bAllOLS'
return matrix Psi =`Psi'
return matrix sPsi =`sPsi'
return matrix ePsi =`ePsi'
return matrix stdvec =`stdvec' // penalized vars only
return scalar lambda =`lambda' // relates to depvar in original units
return scalar slambda =`slambda' // relates to standardized depvar
return scalar lambda0 =`lambda0' // BCH definition of lambda
return scalar rmse =`rmse' // lasso rmse
return scalar r2 =`r2'
return scalar objfn =`objfn'
return scalar c =`c'
return scalar gamma =`gamma'
return scalar gammad =`gammad'
return scalar rmseOLS =`rmseOLS' // post-lasso rmse
return local selected0 `selected0' // all selected/notpen vars INCLUDING NOTPEN (but excl constant)
return local selected `selected' // all selected (penalized) vars EXCLUDING NOTPEN & CONS
return scalar k =`k' // number of all vars in sel/notpen parameter vector INCLUDING CONSTANT
return scalar s0 =`s0' // number of vars selected INCLUDING NOTPEN (but excl constant)
return scalar s =`s' // number of vars selected EXCLUDING NOTPEN & CONS
return scalar p0 =`p0' // number of all vars in original model including constant
return scalar p =`p' // p defined by BCH; excludes notpen & cons
return scalar N_clust =`N_clust'
return scalar N_clust1 =`N_clust1'
return scalar N_clust2 =`N_clust2'
return scalar N =`N'
return scalar center =`center'
// for HAC lasso
return scalar bw =`bw'
return local kernel `kernel'
return local psinegs =`psinegs'
return local psinegvars `psinegvars'
return local clustvar `cluster'
return local robust `robust'
return scalar niter =`niter'
return scalar maxiter =`maxiter'
return scalar npsiiter =`npsiiter'
return scalar maxpsiiter=`maxpsiiter'
return scalar ssnumsim =`ssnumsim'
return scalar supscore =`CCK_ss'
return scalar supscore_p =`CCK_p'
return scalar supscore_cv =`CCK_cv'
return scalar supscore_gamma=`CCK_gamma'
*
end // end _rlasso subroutine
// Subroutine for lassopath
program define _lassopath, rclass sortpreserve
version 13
syntax [anything] , /// anything is actually a varlist but this is more flexible - see below
toest(varlist numeric max=1) /// required `touse' variable (so marksample unneeded)
[ ///
verbose(int 0) ///
consflag(int 1) /// default is constant specified
stdallflag(int 0) /// =1 if lambdas etc. are provided in the standardized metric
dofminus(int 0) /// lost degrees of freedom from FEs
sdofminus(int 0) /// lost degrees of freedom from partialling out ("s"=small; includes the constant)
pminus(int 0) /// omitted/base variables to be excluded from model count
dmflag(int 0) /// data have been demeaned or should be treated as demeaned
notpen_t(string) ///
notpen_o(string) ///
///
/// lambda settings
ALPha(real 1) /// elastic net parameter
Lambda(string) /// overwrite default lambda
lambda2(string) /// overwrite default L2 lambda
LCount(integer 100) ///
LMAX(real 0) ///
LMINRatio(real 1e-4) /// ratio of maximum to minimum lambda
///
TOLOpt(real 1e-10) /// was originally 1e-5 but 1e-10 gives good equiv between xstd and stdloadings
TOLZero(real 1e-4) ///
MAXIter(integer 10000) /// max number of iterations in estimating lasso
fdev(real 1e-5) /// minimum fractional change in deviance for stopping path a la glmnet
devmax(real 0.999) /// maximum fraction of explained deviance for stopping path a la glmnet
///
PLoadings(string) /// set penalty loadings as Stata row vector
ploadings2(string) /// set optional L2 penalty loadings as Stata row vector
///
LGLMnet /// use glmnet scaling for lambda/alpha
sklearn /// use Python sklearn
///
xnames_o(string) /// original names in varXmodel if tempvars used (string so varlist isn't processed)
xnames_t(string) /// temp names in varXmodel
///
sqrt /// square-root lasso
ols ///
///
STDLflag(int 0) /// use standardisation loadings?
stdy(string) ///
stdx(string) ///
///
ADAPTive /// adaptive lasso
ADATheta(real 1) /// gamma paramater for adapLASSO
ADALoadings(string) ///
ADAFirst(string) /// specify first-step ada estimator
///
holdout(varlist) ///
///
NOIC ///
EBICGamma(real -99) /// -99 leads to the default option
///
NODEVCRIT /// do not use deviance as criteria to exit path
]
*** `anything' is actually a varlist but this is more flexible.
// avoids problems with uninitialized vars, unwanted parsing of factor vars, etc.
// ... so immediately rename to `varlist' (usual Stata convention).
local varlist `anything'
*** quietly unless vverbose specified
if `verbose'<2 {
local quietly quietly
}
*
*** count number of obs
_nobs `toest'
local nobs = r(N)
*
*** define various parameters/locals/varlists
local varY `varlist' // now is only element in varlist
local varXmodel `xnames_t' // including notpen if any; may include tempnames
local pen_t : list xnames_t - notpen_t // list of penalized regressors
local pmodel : word count `varXmodel' // # regressors (excl. constant) minus partial
local p0 : word count `varXmodel' // # all regressors incl notpen but excl constant
local p0 = `p0' + `consflag' // ditto but now incl constant
*
*** syntax checks
if (`lminratio'<0) | (`lminratio'>=1) {
di as err "lminratio should be greater than 0 and smaller than 1."
exit 198
}
// txt rather than err in order to stop message unless called by lasso2
if ("`lambda'"!="") & ((`lcount'!=100) | (`lmax'>0) | (`lminratio'!=1e-4)) {
di as txt "lcount | lmax | lminratio option ignored since lambda() is specified."
}
if (`ebicgamma'!=-99) & (`ebicgamma'<0) & (`ebicgamma'>1) {
di as err "ebicgamma must be in the interval [0,1]. Default option is used."
}
*
*** useful flags and counts
local prestdflag =("`stdy'"!="")
local sqrtflag =("`sqrt'"!="")
local olsflag =("`ols'"!="")
local adaflag =("`adaptive'"!="")
local lglmnetflag =("`lglmnet'"!="")
local noicflag =("`noic'"!="")
local nodevcritflag =("`nodevcrit'"!="")
local sklearnflag =("`sklearn'"!="")
if (~`consflag') local noconstant noconstant
*
*** syntax checks
if `sqrtflag' & `lglmnetflag' {
di as err "lglmnet option not supported for square-root lasso"
exit 198
}
if `sklearnflag' {
// check that Stata is at least version 16 and Python is available
cap python query
if _rc==199 {
// Stata 15 or lower
di as err "error - sklearn option requires Stata 16 or higher and a Python-aware Stata"
exit 198
}
else if _rc==111 {
// Stata 16 but no Python
di as err "error - sklearn option requires a Python-aware Stata - see {helpb help python}"
exit 111
}
else if _rc>0 {
// some other error code
di as err "error - sklearn option requires Stata 16 or higher and a Python-aware Stata"
exit 198
}
// sklearn requires lglmnet
if ~`lglmnetflag' {
di as err "error - sklearn option requires lglmnet option"
exit 198
}
// and no custom penalty loadings
if "`ploadings'`ploadings2'"~="" {
di as err "error - sklearn option does not currently support custom penalty loadings"
exit 198
}
}
*** special case - lglmnet with unit loadings (=not prestandardized)
if `lglmnetflag' & "`stdy'"=="" {
tempvar tY
tempname stdy stdx
qui sum `varY' if `toest'
local sdY = r(sd) * sqrt((r(N)-1)/r(N))
local mY = r(mean)
if `consflag' {
qui gen double `tY' = `mY' + (`varY' -`mY') * 1/`sdY' if `toest'
}
else {
qui gen double `tY' = `varY' * 1/`sdY' if `toest'
}
local varY `tY'
mat `stdy' = `sdY'
mat `stdx' = J(1,`p0'-`consflag',1)
local prestd prestd
local prestdflag = 1
local stdlflag = 0
}
*
*********************************************************************************