#this is R version of aenet2 code in Gauss for Zhang Caner paper #generate data first #of iterations ite=10 d=5 #dimesnion of instruments rmse=matrix(0,ite,d) mse=matrix(0,ite,1) bias=matrix(0,ite,d) sel=matrix(0,ite,1) for (jj in (1:ite)){ t=100 co=0.5 rh=co*rep(1,d) rh=matrix(rh,d,1) rh1=sqrt((1-co^2))*rep(1,d) rh1=matrix(rh1,d,1) eta=rnorm(t*d) eta=matrix(eta,t,d) v=rnorm(t*d) v=matrix(v,t,d) e=(eta%*%rh)+(v%*%rh1) #t by 1 e vector omz=diag(1,d,d) kap=0.5 omz[2,1]=kap omz[1,2]=kap z1=chol(omz)%*%matrix(rnorm(t*d),d,t) x=t(z1)#instruments t by d p=rep(1,d^2) p=matrix(p,d,d) z=x%*%p+eta #red form is done theta=c(3,3,3,3,0) theta=matrix(theta,5,1) y=z%*%theta+e #st eqn is done #data generation is done th1hat=solve(t(z)%*%x%*%t(x)%*%z)%*%(t(z)%*%x%*%t(x)%*%y) avg=mean(y-z%*%th1hat) avg1=avg*rep(1,t) shh1=t(((y-z%*%th1hat)-avg1))%*%((y-z%*%th1hat)-avg1)/t #true until here with gauss lat1s=seq(1,1.3,0.1) lat1s=matrix(lat1s,4,1) #putting zeros for matrices to startup. may be we donot need this in R? bico=rep(0,4) bico=matrix(bico,4,1) theta1=rep(0,4*d) theta1=matrix(theta1,d,4) atheta1=rep(0,4*d) atheta1=matrix(atheta1,d,4) atheta2=rep(0,4*d) atheta2=matrix(atheta2,d,4) atheta3=rep(0,4*d) atheta3=matrix(atheta3,d,4) atheta4=rep(0,4*d) atheta4=matrix(atheta4,d,4) atheta=rep(0,4*d) atheta=matrix(atheta,d,4) athetas=rep(0,4*d) athetas=matrix(athetas,d,4) d1=rep(0,4) d2=rep(0,4) d3=rep(0,4) d4=rep(0,4) d5=rep(0,4) d6=rep(0,4) d7=rep(0,4) d1=matrix(d1,4,1) d2=matrix(d2,4,1) d3=matrix(d3,4,1) d4=matrix(d4,4,1) d5=matrix(d5,4,1) d6=matrix(d6,4,1) d7=matrix(d7,4,1) #################################### for (ii in (1:4)){ #change lat1s[1] to lat1s[ii] #(1+(lat1s[1]/t))*((y-z*th)'*x*inv(x'x)*x'*(y-z*th)/shh1)+(lat1s[1]*sumc(abs(th)))+(lat1s[1]*sumc((th.^2)))) #function defn fa=function(th){ upenf=t((y-z%*%th))%*%x%*%solve(t(x)%*%x)%*%t(x)%*%(y-z%*%th)/shh1 penf=lat1s[ii]*sum(abs(th))+lat1s[ii]*sum(th^2) of=(1+(lat1s[ii]/t))*(upenf+penf)} th=c(2.9, 2.9, 3.1, 3.05, 0.1) th=matrix(th,d,1) enetf=optim(th,fa,method="BFGS",hessian=TRUE,control=list(maxit=1)) thetaen=enetf\$par #that 5 by 1 column vector thetaen=matrix(thetaen,d,1) theta1[,ii]=thetaen #setup weights wj=theta1[,ii]^(-4.5) wj=matrix(wj,d,1) #function for aenet, does not allow square bracket in parameter value in function() fanet=function(thetaen){ upenf=t((y-z%*%theta1[,ii]))%*%x%*%solve(t(x)%*%x)%*%t(x)%*%(y-z%*%theta1[,ii])/shh1 penf=lat1s[ii]*sum(wj*abs(theta1[,ii]))+lat1s[ii]*sum(theta1[,ii]^2) of=(1+lat1s[ii]/t)*(upenf+penf)} aentf=optim(theta1[,ii],fanet,method="BFGS",hessian=TRUE,control=list(maxit=1)) atheta1[,ii]=aentf\$par #that 5 by 1 column vector gi=t(x)%*%z w=solve(t(x)%*%x/t) atheta2[,ii]=atheta1[,ii] atheta3[,ii]=atheta1[,ii] f0j=rep(0,d) foj=matrix(f0j,d,1) #I think the outcomes are d by 1, but written as row in output #cannot use rbind to get uneven rows to be binnded, so use matrix to get them in rows for (i in 2:(d-1)) { atheta1[i,ii]=0 yy=atheta1[1:(i-1),ii] yy=matrix(yy,i-1,1) yy1=atheta1[i,ii] yy1=matrix(yy1,1,1) yy2=atheta1[(i+1):d,ii] yy2=matrix(yy2,(d-i),1) atheta[,ii]=rbind(yy,yy1,yy2) sco=t(x)%*%(y-z%*%atheta[,ii]) f0=abs(-t(gi)%*%w%*%sco) f0j[i]=(f0[i]<=lat1s[ii]%*%wj[i]) atheta1[i,ii]=atheta[i,ii]*(1-f0j[i]) if (f0j[i]==1) {atheta3[i,ii]=0} else {atheta3[i,ii]=atheta2[i,ii]} atheta1[i,ii]=atheta2[i,ii] } kk=atheta2[1:(d-1),ii] kk=matrix(kk,d-1,1) athetas[,ii]=rbind(kk,0) scod=t(x)%*%(y-z%*%athetas[,ii]) f0d=abs(t(-gi)%*%w%*%scod) f0jd=(f0d[d]<=lat1s[ii]*wj[d]) atheta3[d,ii]=atheta[d,ii]*(1-f0jd) d3[ii]=(atheta3[4,ii]==0) d4[ii]=(atheta3[5,ii]==0) d5[ii]=(atheta3[3,ii]==0) d6[ii]=(atheta3[2,ii]==0) d7[ii]=(atheta3[1,ii]==0) nz=d-(d3[ii]+d4[ii]+d5[ii]+d6[ii]+d7[ii]);#number of nonzero parameters bico[ii]=t((y-z%*%atheta3[,ii]))%*%x%*%solve(t(x)%*%x)%*%t(x)%*%(y-z%*%atheta3[,ii])/t+((log(t)*nz)/t) } #bic objec. function mbici=which.min(bico)#gets min bic index mlam=lat1s[mbici]#tuning parameter choice sel[jj]=d4[mbici]-d3[mbici]-d5[mbici]-d6[mbici]-d7[mbici]#if nonzero, then shows a mistake, if zero that is correct #create identity matrix, of order d eyed=diag(1,d,d) sig=t(p)%*%omz%*%p+eyed rmse[jj,]=t(((atheta3[,mbici]-theta)^2)) mse[jj]=t((atheta3[,mbici]-theta))%*%sig%*%(atheta3[,mbici]-theta) bias[jj,]=t((atheta3[,mbici]-theta)); #print(bias[jj,]) } print( "Correct percentage:") print(colMeans(sel)) print("MSE") print(colMeans(mse)) print("bias: ") print(colMeans(bias)) print("rmse: ") print(sqrt(colMeans(rmse)))