This markdown file contains attempts to replicate the results in the paper titled “INFERENCE FOR HIGH-DIMENSIONAL SPARSE ECONOMETRIC MODELS” by A. BELLONI, V. CHERNOZHUKOV, AND C. HANSEN.
library(glmnet)## Loading required package: Matrix## Loaded glmnet 4.1-8library(pracma) ##
## Attaching package: 'pracma'## The following objects are masked from 'package:Matrix':
##
## expm, lu, tril, triuset.seed(123)nsim <- 200
n <- 100
p <- 200
alpha <- 1
bx <- 5
ex <- 10
Sigma <- toeplitz(0.5^(0:(p-1)))
C <- chol(Sigma)
beta <- c(1/(1:bx), rep(0, bx), 1/(1:bx), rep(0, p - 3*bx))
eta <- c(1/(1:ex), rep(0, p - ex))
c<- 0.005 # change the C here (do not change more than 0.005)
psi<-10
gamma<- 0.05
#stores alpha coefficient for each simulation
lasso <-post_lasso<- indirect_post<- double_selection<- numeric(nsim)
# Final summary table
df<- data.frame(Estimator=c("Lasso","Post-Lasso","Indirect Post-Lasso","Double selection"),
MeanBias=rep(NA,4),
StdDev=rep(NA,4),
rp=rep(NA,4))# Estimation of SIGMA necessary to find the lambda (based on the defined as the 1-gamma quantile)
# The last argument repe (repetition) has nothing to do with the number of simulations
quantile_part<- function(X,gamma=0.05,repe=100){
dis<- numeric(repe)
#X<- apply(X, 2, function(col) col / norm(col / sqrt(n), type = '2'))
for (l in 1:repe){
#g<- matrix(rnorm(nrow(X)*ncol(X)),nrow=nrow(X))
g<- matrix(rnorm(nrow(X)))
#g<- apply(g, 2, function(col) col / norm(col / sqrt(n), type = '2')) # normalize g as well from MC_TE_simulateLamdasqrtLASSO.m
dis[l]<- max(abs(t(X)%*%g)) # distribution
#dis[l]<- nrow(X)*max(abs(t(X)%*%g)) # distribution
}
#cat("Quantile value is ",quantile(dis,probs=1-gamma), '\n' )
return(as.numeric(quantile(dis,probs=1-gamma)))
}
## Generate data (X) based on the normalization
generateX<-function(n,p,C){
X<-matrix(rnorm(n*p), ncol = p)%*%C # generate x
# normalize each column of x
for (k in 1:ncol(X)){
X[,k]<- X[,k] / norm(X[,k]/sqrt(n),type='2') # page 3 condition asm and comment 4.1
#X[,k]<- X[,k] / norm(X[,k],type='2') # page 3 condition asm and comment 4.1
}
#X <- apply(X, 2, function(col) col / norm(col / sqrt(n), type = '2'))
return(X)
}
# ols for sigma estimation (from LassoShooting.m)
ols_sigma<- function(X,y){
beta_esti<-solve((t(X) %*% X + (0.01* eye(p))))%*% t(X) %*% y
return(sd(y- X%*% beta_esti)) # residual sd
#return(29)
}
ols_esti<- function(X,y){
if(!is.matrix(X)){ # if one column is selected
X=matrix(X,ncol = 1)
}
post_lasso_model <- lm(y ~ X) #with intercept
coef_post_lasso <- coef(summary(post_lasso_model))
return(coef_post_lasso[2,1])
}
#calculates the optimal sigma, as defined in Appendix A of Belloni (2011)
#K, the upper bound on the number of iterations is defined as >1, hence is set at 100
#c and gamma as well are given generic values of 1.1 and 0.05 respectively
# v (tolerance) is set to 1e-3 (we use this rather than number of iterations K, though K is kept as an argument)
estimate_lambda<- function(X,y,c=1.1,psi=1,gamma=0.05,K=100,v=1e-3){
#sigma I0 defined, I0 is a set of regressors that is included in the model (from the beginning of Appendix A)
sigma<-array(NA,dim=K)
# initial OLS
sigma[1]<- psi *ols_sigma(X,y)
k=2
cond=TRUE
qpart<- quantile_part(X,gamma)
#cat('Quantile is ', qpart , '\n')
while(cond){
#find the lambda for lasso
lambda_est<- 2* c * sigma[k-1] * qpart
#cat('Lambda at ', k, ' is ', lambda_est , '\n')
#fit the lasso
lasso_reg <- glmnet(X, y,lambda=lambda_est,intercept = F ) #change lambda to x dependent lambda when it works
#cat(dim(cbind(d,X)))
#cat(length(as.numeric(coef(lasso_reg)[-1])!=0) , '\n ')
lasso_esti<-X %*% as.numeric(coef(lasso_reg)[-1])
lasso_resi<- y-lasso_esti
sigma[k]<- sd(lasso_resi) #sigma^2 is Q(B), where Q(B) is variance of the residuals
#stop the iteration when sigma starts to converge or when the number of iterations exceeded
if(abs(sigma[k-1]-sigma[k]) < v | (k > K)){
cond=FALSE
#cat('Sigma is ' , sigma[k] , '\n')
#cat('Lambda is ' , lambda_est , '\n')
#print(sigma[!is.na(sigma)])
}else{
k=k+1
}
}
return(lambda_est)
}
#POST Lamda estimation
estimate_lambda_post_lasso<- function(X,y,c=1.1,psi=1,gamma=0.05,K=100,v=1e-3){
sigma<-array(NA,dim=K)
# initial OLS
sigma[1]<- psi *ols_sigma(X,y)
k=2
cond=TRUE
qpart<- quantile_part(X,gamma)
while(cond){
# Find the lambda for lasso
lambda_est<- 2* c * sigma[k-1] * qpart
#cat('Lambda at ', k, ' is ', lambda_est , '\n')
#cat(dim(X),length(y))
# Fit the lasso
lasso_reg <- glmnet(X, y,lambda=lambda_est,intercept = F) # change lambda to x dependent lambda when it works
lasso_coef<- coef(lasso_reg)[-1]
lasso_esti <-X %*% as.numeric(coef(lasso_reg)[-1])
non_zero_coeffs<- which(coef(lasso_reg)[-1]!=0)
S<-length(non_zero_coeffs)
#non_zero_coeffs<- non_zero_coeffs[non_zero_coeffs!=2]
if(S>nrow(X)){
sigma[k]<-ols_sigma(X[,non_zero_coeffs],y)
}else{
fit<-lm(y~X[,non_zero_coeffs] )
var_resi<- var(y-fit$fitted.values)*n/(n-S)
sigma[k]<- sqrt(var_resi)
}
cat('Sigma', sigma[k],'\n')
if(abs(sigma[k-1]-sigma[k]) < v){
cond=FALSE
}else{
k=k+1
}
}
#cat(sigma[1:10])
return(lambda_est)
}#Monte Carlo Simulation!
for (i in 1:nsim) {
set.seed(i)
X <- generateX(n,p,C)
d <- X%*%eta + rnorm(n)
X_d <- cbind(d, X) #this combines the X and d, such that both alpha and beta are estimated
y <- alpha*d + X%*%beta +rnorm(n)
#first, find the optimal sigma that can be used for the lambda
esti_lambda <-estimate_lambda(X=X,y=y,c=c,psi = psi)
#cat("Estimated lambda",esti_lambda , '\n')
## 1. Iterative LASSO
#lasso with X_d as X such that X and d are both included, penalty factor such that only beta is penalized, not alpha
lasso_reg <- glmnet(X_d, y,
lambda=esti_lambda,
#lambda=0.6,
intercept = F,
penalty.factor = c(0,rep(1,p)))
coefficients_lasso <- as.vector(coef(lasso_reg))
lasso[i] <- coefficients_lasso[2] #stores the alpha coefficient (as the first coefficient denotes the intercept)
#selected_lasso_reg <- which(coefficients_lasso_reg != 0) #exclude intercept step one
## 2. LASSO-POST
selected_indices <- which(coefficients_lasso[-1] != 0)
X_selected <- X_d[, selected_indices]
post_lasso[i] <- ols_esti(X_selected,y)
## 3.Indirect post
#lasso with X_d as X such that X and d are both included, penalty factor such that only beta is penalized, not alpha
lasso_indirect <- glmnet(X, d,lambda=esti_lambda,intercept = F)
coefficients_lasso_indirect <- as.vector(coef(lasso_indirect ))
non_zero_locations <- which(coefficients_lasso_indirect[-1]!=0) #stores the alpha coefficient (as the first coefficient denotes the intercept)
X_selected <- X_d[, non_zero_locations]
indirect_post[i] <- ols_esti(X_selected,y)
## 4. Double selection
# Step 1 first, find the optimal sigma for y on X
esti_lambda <-estimate_lambda(X=X,y=y,c=c,psi = psi)
lasso_reg <- glmnet(X, y,lambda=esti_lambda,intercept = F)
step1 <- as.vector(coef(lasso_reg))
non_zero_locations_step1 <- which(step1[-1]!=0) #stores the alpha coefficient (as the first coefficient denotes the intercept)
# Step 2 first, find the optimal sigma for y on X
esti_lambda<- estimate_lambda(X,d,c=c)
lasso_reg <- glmnet(X, d,lambda=esti_lambda,intercept = F)
step2 <- as.vector(coef(lasso_reg))
non_zero_locations_step2 <- which(step2[-1]!=0)
# union of step 1 and 2 coeffieicnts
unique_support<- unique(c(non_zero_locations_step1,non_zero_locations_step2))
double_selection[i]<- ols_esti(cbind(d,X[,unique_support]),y)
cat('Finished ', i , ' iterations \n')
if (i==nsim){
# lasso
df[1,2]<- mean(lasso -alpha)
df[1,3]<- sd(lasso)
Zlasso<-(lasso-alpha)/sd(lasso)
df[1,4]<- sum(abs(Zlasso)>1.96)/nsim
#post-lasso
df[2,2]<- mean(post_lasso -alpha)
df[2,3]<- sd(post_lasso)
Zlasso_post<-(post_lasso-alpha)/sd(post_lasso)
df[2,4]<- sum(abs(Zlasso_post)>1.96)/nsim
#indirect post-lasso
df[3,2]<- mean(indirect_post -alpha)
df[3,3]<- sd(indirect_post)
Zlasso_indirect_post<-(indirect_post-alpha)/sd(indirect_post)
df[3,4]<- sum(abs(Zlasso_indirect_post)>1.96)/nsim
# double selection
df[4,2]<- mean(double_selection -alpha)
df[4,3]<- sd(double_selection)
Zlasso_double_selection<-(double_selection-alpha)/sd(double_selection)
df[4,4]<- sum(abs(Zlasso_double_selection)>1.96)/nsim
df[,2:4]<- round(df[2:4],digits = 4)
print(df)
}
} ## Finished 1 iterations
## Finished 2 iterations
## Finished 3 iterations
## Finished 4 iterations
## Finished 5 iterations
## Finished 6 iterations
## Finished 7 iterations
## Finished 8 iterations
## Finished 9 iterations
## Finished 10 iterations
## Finished 11 iterations
## Finished 12 iterations
## Finished 13 iterations
## Finished 14 iterations
## Finished 15 iterations
## Finished 16 iterations
## Finished 17 iterations
## Finished 18 iterations
## Finished 19 iterations
## Finished 20 iterations
## Finished 21 iterations
## Finished 22 iterations
## Finished 23 iterations
## Finished 24 iterations
## Finished 25 iterations
## Finished 26 iterations
## Finished 27 iterations
## Finished 28 iterations
## Finished 29 iterations
## Finished 30 iterations
## Finished 31 iterations
## Finished 32 iterations
## Finished 33 iterations
## Finished 34 iterations
## Finished 35 iterations
## Finished 36 iterations
## Finished 37 iterations
## Finished 38 iterations
## Finished 39 iterations
## Finished 40 iterations
## Finished 41 iterations
## Finished 42 iterations
## Finished 43 iterations
## Finished 44 iterations
## Finished 45 iterations
## Finished 46 iterations
## Finished 47 iterations
## Finished 48 iterations
## Finished 49 iterations
## Finished 50 iterations
## Finished 51 iterations
## Finished 52 iterations
## Finished 53 iterations
## Finished 54 iterations
## Finished 55 iterations
## Finished 56 iterations
## Finished 57 iterations
## Finished 58 iterations
## Finished 59 iterations
## Finished 60 iterations
## Finished 61 iterations
## Finished 62 iterations
## Finished 63 iterations
## Finished 64 iterations
## Finished 65 iterations
## Finished 66 iterations
## Finished 67 iterations
## Finished 68 iterations
## Finished 69 iterations
## Finished 70 iterations
## Finished 71 iterations
## Finished 72 iterations
## Finished 73 iterations
## Finished 74 iterations
## Finished 75 iterations
## Finished 76 iterations
## Finished 77 iterations
## Finished 78 iterations
## Finished 79 iterations
## Finished 80 iterations
## Finished 81 iterations
## Finished 82 iterations
## Finished 83 iterations
## Finished 84 iterations
## Finished 85 iterations
## Finished 86 iterations
## Finished 87 iterations
## Finished 88 iterations
## Finished 89 iterations
## Finished 90 iterations
## Finished 91 iterations
## Finished 92 iterations
## Finished 93 iterations
## Finished 94 iterations
## Finished 95 iterations
## Finished 96 iterations
## Finished 97 iterations
## Finished 98 iterations
## Finished 99 iterations
## Finished 100 iterations
## Finished 101 iterations
## Finished 102 iterations
## Finished 103 iterations
## Finished 104 iterations
## Finished 105 iterations
## Finished 106 iterations
## Finished 107 iterations
## Finished 108 iterations
## Finished 109 iterations
## Finished 110 iterations
## Finished 111 iterations
## Finished 112 iterations
## Finished 113 iterations
## Finished 114 iterations
## Finished 115 iterations
## Finished 116 iterations
## Finished 117 iterations
## Finished 118 iterations
## Finished 119 iterations
## Finished 120 iterations
## Finished 121 iterations
## Finished 122 iterations
## Finished 123 iterations
## Finished 124 iterations
## Finished 125 iterations
## Finished 126 iterations
## Finished 127 iterations
## Finished 128 iterations
## Finished 129 iterations
## Finished 130 iterations
## Finished 131 iterations
## Finished 132 iterations
## Finished 133 iterations
## Finished 134 iterations
## Finished 135 iterations
## Finished 136 iterations
## Finished 137 iterations
## Finished 138 iterations
## Finished 139 iterations
## Finished 140 iterations
## Finished 141 iterations
## Finished 142 iterations
## Finished 143 iterations
## Finished 144 iterations
## Finished 145 iterations
## Finished 146 iterations
## Finished 147 iterations
## Finished 148 iterations
## Finished 149 iterations
## Finished 150 iterations
## Finished 151 iterations
## Finished 152 iterations
## Finished 153 iterations
## Finished 154 iterations
## Finished 155 iterations
## Finished 156 iterations
## Finished 157 iterations
## Finished 158 iterations
## Finished 159 iterations
## Finished 160 iterations
## Finished 161 iterations
## Finished 162 iterations
## Finished 163 iterations
## Finished 164 iterations
## Finished 165 iterations
## Finished 166 iterations
## Finished 167 iterations
## Finished 168 iterations
## Finished 169 iterations
## Finished 170 iterations
## Finished 171 iterations
## Finished 172 iterations
## Finished 173 iterations
## Finished 174 iterations
## Finished 175 iterations
## Finished 176 iterations
## Finished 177 iterations
## Finished 178 iterations
## Finished 179 iterations
## Finished 180 iterations
## Finished 181 iterations
## Finished 182 iterations
## Finished 183 iterations
## Finished 184 iterations
## Finished 185 iterations
## Finished 186 iterations
## Finished 187 iterations
## Finished 188 iterations
## Finished 189 iterations
## Finished 190 iterations
## Finished 191 iterations
## Finished 192 iterations
## Finished 193 iterations
## Finished 194 iterations
## Finished 195 iterations
## Finished 196 iterations
## Finished 197 iterations
## Finished 198 iterations
## Finished 199 iterations
## Finished 200 iterations
## Estimator MeanBias StdDev rp
## 1 Lasso 0.7053 0.0886 1.000
## 2 Post-Lasso 0.6819 0.0759 1.000
## 3 Indirect Post-Lasso 0.3549 0.2024 0.435
## 4 Double selection -0.0179 0.1148 0.055