# solved of esimation equation, with R function fsolve(pracma) # iterated, fewer number of equations #library(pracma) library(neldermead) #'iter' qr_gmm <- function(X, Y, tau=0.5, maxiter = 1e2,ftol=1e-2,asy=0,method=2) { opt2 <- optimset(TolX=ftol,TolFun=ftol,MaxIter=100,MaxFunEvals=300,Display='off') t1 = proc.time() X <- as.matrix(X) Y <- as.matrix(Y) n <- nrow(Y) r <- ncol(Y) p <- ncol(X) mX = colMeans(X) mY = colMeans(Y) Yc <- as.matrix(scale(Y, center = T, scale = FALSE)) Xc <- as.matrix(scale(X, center = T, scale = FALSE)) W = cbind(rep(1,n),X) sigX <- cov(Xc) * (n-1)/n ModelOutput=list() m1= rq(Y~X,tau = tau) alpha <- m1$coefficients[1] beta = matrix(m1$coefficients[-1],p,1) muX = mX SigmaX = sigX theta0 = c(alpha,muX,beta,SigmaX) # #========= h<-function(theta){ alpha=theta[1] numpar=1 muX=theta[1:p+numpar] numpar=numpar+p beta = matrix(theta[1:p+numpar],p,1) numpar = numpar+p SigmaX = inv.vech(theta[1:(p*(p+1)/2)+numpar]) ones=matrix(1,n,1) SX <- t(X-ones%*%muX)%*%(X-ones%*%muX)/n h1 = t(W)%*%((Y<(X%*%beta+alpha))-tau)/n h2 = muX-mX h3 = vech( SigmaX-SX) c(h1,h2,h3) } #====== f<-function(theta) sum(h(theta)^2) #m1=fminsearch(f=f,x0=theta0) if(method==1) {m0=nlm(f=f,p=theta0,gradtol = ftol,iterlim=maxiter,print.level = 1);theta=m0$estimate} else if(method==2) {m0=neldermead::fminsearch(fun=f,x0=theta0,opt2);theta=(m0$optbase)$xopt} else if(method==3) {m0=optim(par = theta0,fn = f,method = "Nelder-Mead",control = list(fnscale=ftol,maxit=maxiter));theta=m0$par} else {m0=pracma::fminsearch(f=f,x0=theta0,maxiter=maxiter,tol=ftol);theta=m0$xval} #f(m0$estimate) #f(m1$xval) #theta=m0$xval alpha<-theta[1] numpar=1 muX=theta[1:p+numpar] numpar=numpar+p beta = matrix(theta[1:p+numpar],p,1) numpar = numpar+p SigmaX = inv.vech(theta[1:(p*(p+1)/2)+numpar]) DeltaM=0 for(i in 1:n){ g1 = t(W[i,])*(as.numeric(Y[i,]<(X[i,]%*%beta+alpha))-tau) g2 = muX-X[i,] g3 = vech(SigmaX-(X[i,]-muX)%*%t(X[i,]-muX)) gi = c(g1,g2,g3) DeltaM = DeltaM+gi%*%t(gi) } DeltaM=DeltaM/n #print(diag(DeltaM)) DeltaM_inv=chol2inv(sechol(DeltaM+diag(1e-3,ncol(DeltaM)))) #print(eigen(DeltaM_inv,only.values = TRUE)) #print(diag(DeltaM_inv)) fnew<-function(theta) {hval=h(theta); hval%*%DeltaM_inv%*%hval} # #m1=fminsearch(f=fnew,x0=theta,maxiter=100,tol=10-3,dfree = FALSE) if(method==1) {m1=nlm(f=fnew,p=theta,gradtol = ftol,iterlim=maxiter);theta=m1$estimate} else if(method==2){m1=neldermead::fminsearch(fun=fnew,x0=theta,opt2);theta=(m1$optbase)$xopt; logobj=(m1$optbase)$fopt} else if(method==3) {m1=optim(par = theta,fn = fnew,method = "Nelder-Mead",control = list(fnscale=ftol,maxit=maxiter));theta=m1$par} else {m1=pracma::fminsearch(f=fnew,x0=theta,maxiter=maxiter,tol=ftol);theta=m1$xval} alpha<-theta[1] numpar=1 muX=theta[1:p+numpar] numpar=numpar+p beta = matrix(theta[1:p+numpar],p,1) numpar = numpar+p SigmaX = inv.vech(theta[1:(p*(p+1)/2)+numpar]) numpar=r + p * r + p + p * (p + 1) / 2; #Z = Y-(X%*%beta+alpha) #logobj = log(sum(Z*(tau-(Z<0)))/n) BIC = logobj+log(n)*numpar/(2*n) ModelOutput$alpha = alpha ModelOutput$beta = beta #ModelOutput$Gamma = Gamma ModelOutput$paramNum = r + p * r +p + p * (p + 1) / 2; #ModelOutput$Gamma0 = Gamma0 #ModelOutput$eta = eta #ModelOutput$SigX = SigX #ModelOutput$iter = iter #ModelOutput$paramNum = r + u * r +p + p * (p + 1) / 2; #ModelOutput$Gamma0 = Gamma0 #ModelOutput$eta = eta #ModelOutput$SigX = SigX #ModelOutput$iter = iter #ModelOutput$Omega = Omega #ModelOutput$Omega0 = Omega0 #ModelOutput$sigYcX = sigYcX ModelOutput$logobj = logobj ModelOutput$BIC = BIC ModelOutput$paramNum = r + p* r +p + p * (p + 1) / 2 ModelOutput$n = n class(ModelOutput) <- "qev_gmm" ModelOutput$fit_time=(proc.time()-t1)[3] return(ModelOutput) }