################################################################################
#---------------------@Taining Wang Simulation Code----------------------------#
#-----------Local within estimation on varying coefficient SF panel model------#
#-----------------------With Feng Yao and Jun Cai------------------------------#
#------------ISEM, Capital University of Economics and Business----------------#
#----------------------First version: Oct 13, 2023-----------------------------#
#------------------------Last Update: Jan 15, 2026-----------------------------#
################################################################################

library(minqa)
library(compiler)

#---[1] Within transformation over t for a matrix x in panel data
within.tran<-function(x,n,t,constrain=FALSE){
  
  nT<-n*t
  x<-as.matrix(x)
  p<-ncol(x)
  Md.x<-matrix(0,nT,p)
  if(constrain==FALSE){
    for (j in 1:p) {
      xj.nt<-matrix(x[,j],n,t,byrow=T)
      xj.n.bar<-rowMeans(xj.nt)
      xj.n.diff<-xj.nt-xj.n.bar
      Md.x[,j]<-matrix(t(xj.n.diff),nT,1)
    }
  }else{
    for (j in 1:p) {
      xj.nt<-matrix(x[,j],n,t,byrow=T)
      xj.n.bar<-rowSums(xj.nt)
      xj.n.sum<-rep(0,n)
      for (jj in 1:n) {
        xj.n.sum[jj]<-(n-1)*xj.n.bar[jj]-sum(xj.n.bar[-jj]) 
      }
      
      xj.n.diff<-xj.nt-(xj.n.sum/nT)
      Md.x[,j]<-matrix(t(xj.n.diff),nT,1)
    }    
  }
  return(Md.x)
}

#---[2] LCLS model to estimate V(eps|x) in inference step
lcls<-cmpfun(function(y,x,x.eval,C.h=NULL,lono='no',ckertype='gaussian'){
  y<-as.matrix(y)
  x<-as.matrix(x)
  x.eval<-as.matrix(x.eval)
  n.eval <- nrow(x.eval)
  n <- length(y)
  q <- ncol(x)
  q1 <- q+1
  ones<-rep(1,n)
  
  if(ckertype=="gaussian") {
    kk <- function(g) (1/sqrt(2*pi))*exp(-0.5*g^2)
    Rk<-1/(2*sqrt(pi))
    if(is.null(C.h)==TRUE){
      h<-1.059*apply(x,2,sd)*n^(-1/(4+q))
    }else{
      h<-C.h*apply(x,2,sd)*n^(-1/(4+q))
    }
  }else{
    if(ckertype=='uniform'){
      kk<-function(g) ifelse(abs(g)<=1,0.5,0)
      Rk<-1/2
      if(is.null(C.h)==TRUE){
        h<-1.843*apply(x,2,sd)*n^(-1/(4+q))
      }else{
        h<-C.h*apply(x,2,sd)*n^(-1/(4+q))
      }
    }else{
      if(ckertype=="epanechnikov") {
        kk <- function(g) ifelse(abs(g)<=1, 0.75*(1-g^2), 0)
        Rk<-3/5
        if(is.null(C.h)==TRUE){
          h<-2.345*apply(x,2,sd)*n^(-1/(4+q))
        }else{
          h<-C.h*apply(x,2,sd)*n^(-1/(4+q))
        }
      }else{
        if(ckertype=="biweight") {
          kk <- function(g) ifelse(abs(g)<=1, (15/16)*(1-g^2)^2, 0)
          Rk<-5/7
          if(is.null(C.h)==TRUE){
            h<-2.778*apply(x,2,sd)*n^(-1/(4+q))
          }else{
            h<-C.h*apply(x,2,sd)*n^(-1/(4+q))
          }
        }else{
          if(ckertype=="triweight") {
            kk <- function(g) ifelse(abs(g)<=1, (35/32)*(1-g^2)^3, 0)
            Rk<-350/429
            if(is.null(C.h)==TRUE){
              h<-3.154*apply(x,2,sd)*n^(-1/(4+q))
            }else{
              h<-C.h*apply(x,2,sd)*n^(-1/(4+q))
            }
          }
        }
      }
    }
  }
  
  thetahat <- matrix(0,n.eval,q1)
  for(j in 1:n.eval){ 
    
    K<-ones
    for (i in 1:q){
      v <- (x[,i] - x.eval[j,i])/h[i]
      K<-K*kk(v)
    }
    if(lono=='yes') K[i]=0
    thetahat[j,1] <- sum(y*K)/sum(K) 
    
  }
  
  result.list<-list(fitted.y=thetahat[,1],h=h,Rk=Rk)
  return(result.list)
})

#---[3] LLLS for pool model
llls<- function(y,x,x.eval,C.h=NULL,lono='no',ckertype='gaussian'){
  
  x<-as.matrix(x)
  x.eval<-as.matrix(x.eval)
  n.eval<-nrow(x.eval)
  y<-as.matrix(y)
  n <- length(y)
  q <- ncol(x)
  q1 <- q+1
  ones <- rep(1,n)
  
  if(ckertype=="gaussian") {
    kk <- function(g) (1/sqrt(2*pi))*exp(-0.5*g^2)
    if(is.null(C.h)==TRUE){
      h<-1.059*apply(x,2,sd)*n^(-1/(4+q))
    }else{
      h<-C.h*apply(x,2,sd)*n^(-1/(4+q))
    }
  }else{
    if(ckertype=='uniform'){
      kk<-function(g) ifelse(abs(g)<=1,0.5,0)
      if(is.null(C.h)==TRUE){
        h<-1.843*apply(x,2,sd)*n^(-1/(4+q))
      }else{
        h<-C.h*apply(x,2,sd)*n^(-1/(4+q))
      }
    }else{
      if(ckertype=="epanechnikov") {
        kk <- function(g) ifelse(abs(g)<=1, 0.75*(1-g^2), 0)
        if(is.null(C.h)==TRUE){
          h<-2.345*apply(x,2,sd)*n^(-1/(4+q))
        }else{
          h<-C.h*apply(x,2,sd)*n^(-1/(4+q))
        }
      }else{
        if(ckertype=="biweight") {
          kk <- function(g) ifelse(abs(g)<=1, (15/16)*(1-g^2)^2, 0)
          if(is.null(C.h)==TRUE){
            h<-2.778*apply(x,2,sd)*n^(-1/(4+q))
          }else{
            h<-C.h*apply(x,2,sd)*n^(-1/(4+q))
          }
        }else{
          if(ckertype=="triweight") {
            kk <- function(g) ifelse(abs(g)<=1, (35/32)*(1-g^2)^3, 0)
            if(is.null(C.h)==TRUE){
              h<-3.154*apply(x,2,sd)*n^(-1/(4+q))
            }else{
              h<-C.h*apply(x,2,sd)*n^(-1/(4+q))
            }
          }
        }
      }
    }
  }
  thetahat <- matrix(0,n.eval,q1)
  for(i in 1:n.eval){ 
    
    X<- matrix(1,n,q+1)
    K <- ones
    for (ii in 1:q){
      v <- (x[,ii] - x.eval[i,ii])/h[ii]
      v.x <- (x[,ii] - x.eval[i,ii])
      K <- K*kk(v)
      X[,ii+1] <- v.x
    }
    if(lono=='yes') K[i]=0
    KX<-K*X
    XK<-t(KX)
    thetahat[i,]<-t((chol2inv(chol(XK%*%X)))%*%(XK%*%y))
  }
  
  result.list<-list(fitted.y=thetahat[,1],gradient.x=thetahat[,2:q1],h=h)
  return(result.list)
}

#---[3] kernel density estimator
npuden<-cmpfun(function(x,x.eval=NULL,h=NULL,lono='no',ckertype="gaussian"){
  
  
  x <- as.matrix(x)
  q <- ncol(x)
  n <- nrow(x)
  ones<-rep(1,n)
  if(is.null(x.eval)) x.eval=x else x.eval<-as.matrix(x.eval)
  n.eval<-nrow(x.eval)
  
  
  if(ckertype=="gaussian") {
    kk <- function(g) (1/sqrt(2*pi))*exp(-0.5*g^2)
    if(is.null(h)==TRUE){
      h<-1.059*apply(x,2,sd)*n^(-1/(4+q))
    }else{
      h<-h
    }
  }else{
    if(ckertype=='uniform'){
      kk<-function(g) ifelse(abs(g)<=1,0.5,0)
      if(is.null(h)==TRUE){
        h<-1.843*apply(x,2,sd)*n^(-1/(4+q))
      }else{
        h<-h
      }
    }else{
      if(ckertype=="epanechnikov") {
        kk <- function(g) ifelse(abs(g)<=1, 0.75*(1-g^2), 0)
        if(is.null(h)==TRUE){
          h<-2.345*apply(x,2,sd)*n^(-1/(4+q))
        }else{
          h<-h
        }
      }else{
        if(ckertype=="biweight") {
          kk <- function(g) ifelse(abs(g)<=1, (15/16)*(1-g^2)^2, 0)
          if(is.null(h)==TRUE){
            h<-2.778*apply(x,2,sd)*n^(-1/(4+q))
          }else{
            h<-h
          }
        }else{
          if(ckertype=="triweight") {
            kk <- function(g) ifelse(abs(g)<=1, (35/32)*(1-g^2)^3, 0)
            if(is.null(h)==TRUE){
              h<-3.154*apply(x,2,sd)*n^(-1/(4+q))
            }else{
              h<-h
            }
          }
        }
      }
    }
  }
  
  
  h.prod <- prod(h)
  K.matrix <- matrix(0,n,n.eval) 
  #--- n.eval by n.eval kernel matrix
  for (j in 1:n.eval){
    
    KK <- ones
    for (i in 1:q){
      v <- (x[,i] - x.eval[j,i])/h[i]
      KK<- KK*kk(v)
    }
    if(lono=='yes') KK[j]=0
    K.matrix[,j]<-KK
  }
  fhat<-(n*h.prod)^(-1)*colSums(K.matrix)
  
  result.list <- list(fhat=fhat,h=h,K.matrix=K.matrix,ckertype=ckertype)
  return(result.list)
})

#---[4.1] DGP: with intercept function \phi(z); FE normalization can be imposed by
#   setting fe.normalize=T, or not imposed by setting fe.normalize=F
DGP<-cmpfun(function(dgp.type=1,grid.length=50,n,t,c0=1,sigma.v2=1,a=-2,b=2,p=2,
                     sigma.fe=0.5,fe.normalize=FALSE){
  
  nT<-n*t
  
  x1.nt<-matrix(runif(nT,a,b),n,t)
  x2.nt<-matrix(runif(nT,a,b),n,t)
  x3.nt<-matrix(runif(nT,a,b),n,t)
  
  z.nt<-matrix(runif(nT,a,b),n,t)
  
  x1<-as.vector(t(x1.nt))
  x2<-as.vector(t(x2.nt))
  x3<-as.vector(t(x3.nt))
  
  z<-as.vector(t(z.nt))
  
  z.eval<-seq(a,b,length=grid.length)
  
  if(fe.normalize){
    alpha_1<-(c0*(rowMeans(x1.nt+z.nt))+rnorm(n,0,sigma.fe))[-1]
    alpha<-c(-sum(alpha_1),alpha_1)
  }else{
    alpha<-(c0*(rowMeans(x1.nt+z.nt))+rnorm(n,0,sigma.fe))
  }
  
  Dalpha<-rep(alpha,each=t)
  
  if(dgp.type==1){
    
    b0.z<-0.5*z^2; b0.z.eval<-0.5*z.eval^2;
    b1.z<-1+z-z^2; b1.z.eval<-1+z.eval-z.eval^2; 
    b2.z<-sin(0.5*pi*z); b2.z.eval<-sin(0.5*pi*z.eval);
    b3.z<-z^3; b3.z.eval<-z.eval^3;
    
    b0.z.der.eval<-z.eval
    b1.z.der.eval<-1-2*z.eval
    b2.z.der.eval<-cos(0.5*pi*z.eval)*0.5*pi
    b3.z.der.eval<-3*z.eval^2
    
  }else{
    if(dgp.type==2){
      
      b0.z<-cos(z)-1; b0.z.eval<-cos(z.eval)-1;
      b1.z<-exp(z)+1; b1.z.eval<-exp(z.eval)+1; 
      b2.z<-z^3; b2.z.eval<-z.eval^3;
      b3.z<-z^2; b3.z.eval<-z.eval^2;
      
      b0.z.der.eval<--sin(z.eval)
      b1.z.der.eval<-exp(z.eval)
      b2.z.der.eval<-3*z.eval^2
      b3.z.der.eval<-2*z.eval
    }
  }
  
  #---generating error
  v<-rnorm(nT,0,sqrt(sigma.v2))
  
  if(p==1){
    
    y<-b0.z+x1*b1.z+Dalpha+v
    x<-cbind(x1)
    
    fun<-cbind(b0.z,b1.z)
    fun.eval<-cbind(b0.z.eval,b1.z.eval)
    fun.der.eval<-cbind(b0.z.der.eval,b1.z.der.eval)
    
    
  }else{
    if(p==2){
      
      y<-b0.z+x1*b1.z+x2*b2.z+Dalpha+v
      x<-cbind(x1,x2)
      
      fun<-cbind(b0.z,b1.z,b2.z)
      fun.eval<-cbind(b0.z.eval,b1.z.eval,b2.z.eval)
      fun.der.eval<-cbind(b0.z.der.eval,b1.z.der.eval,b2.z.der.eval)
    }else{
      if(p==3){
        y<-b0.z+x1*b1.z+x2*b2.z+x3*b3.z+Dalpha+v
        x<-cbind(x1,x2,x3)
        
        fun<-cbind(b0.z,b1.z,b2.z,b3.z)
        fun.eval<-cbind(b0.z.eval,b1.z.eval,b2.z.eval,b3.z.eval)
        fun.der.eval<-cbind(b0.z.der.eval,b1.z.der.eval,b2.z.der.eval,b3.z.der.eval)
      }
    }
  }
  
  
  return.list<-list(y=y,x=x,z=z,z.eval=z.eval,grid.length=grid.length,c0=c0,
                    fun=fun,fun.eval=fun.eval,fun.der.eval=fun.der.eval,n=n,t=t,nT=nT,
                    fe.normalize=fe.normalize)
  
  return(return.list)    
})

#---[5.1] local-within local linear estimator for one-way VCM: X'b(Z)
#   Note: it returns derivative of b0(z) & function and derivative of b0(z),...,bp(z)
vcm.lwll.step1<-cmpfun(function(y,x,z,z.eval=NULL,n,t,h,ckertype='gaussian',lono=F,Print=T){
  
  y= as.matrix(y); x= as.matrix(x);  z= as.matrix(z);  nT=n*t; h<-as.numeric(h)
  
  if(is.null(z.eval)) z.eval=z else z.eval=as.matrix(z.eval)
  nT.eval<-nrow(z.eval)
  p<- ncol(x)
  q<- ncol(z)
  
  ones<-rep(1,nT)
  ones.T<-rep(1,t)
  
  X<-cbind(ones,x); p1<-p+1
  
  if(ckertype=="gaussian") {
    k <- function(g) (1/sqrt(2*pi))*exp(-0.5*g^2)
    Rk<-1/(2*sqrt(pi))
  }else{
    if(ckertype=='uniform'){
      k<-function(g) ifelse(abs(g)<=1,0.5,0)
      Rk<-1/2
    }else{
      if(ckertype=="epanechnikov") {
        k <- function(g) ifelse(abs(g)<=1, 0.75*(1-g^2), 0)
        Rk<-3/5
      }else{
        if(ckertype=="biweight") {
          k <- function(g) ifelse(abs(g)<=1, (15/16)*(1-g^2)^2, 0)
          Rk<-5/7
        }else{
          if(ckertype=="triweight") {
            k <- function(g) ifelse(abs(g)<=1, (35/32)*(1-g^2)^3, 0)
            Rk<-350/429
          }
        }
      }
    }
  }
  
  #---[1.1] evaluating at grid point
  theta.hat<-matrix(0,nT.eval,p)
  theta.der.hat<-matrix(0,nT.eval,p1*q)
  for(i in 1:nT.eval){ 
    
    if(Print) cat(paste("\rStep 1 Estimation at",i,'of',nT.eval))
    R<-x
    K<-ones
    for (jj in 1:q){
      dz<-(z[,jj]-z.eval[i,jj])
      v.z<-dz/h[jj]
      K<-K*k(v.z)
      
      for (jjj in 1:p1){
        x.dz <- X[,jjj]*dz
        R <- cbind(R,x.dz)
      }
    }
    
    if(lono) K[i]<-0
    
    #---[a]: estimating b(z)
    K.Tn<-matrix(K,t,n)
    
    K.bar.n<-colSums(K.Tn)
    K.bar.nT<-kronecker(K.bar.n,ones.T)
    w.x.a<-K/K.bar.nT
    w.x.a.Tn<-matrix(w.x.a,t,n)
    
    #---[1] transform Y
    Y.Tn<-matrix(y,t,n)
    y.bar.a<-kronecker(colSums(Y.Tn*w.x.a.Tn),ones.T)
    
    Y.til<-y-as.matrix((y.bar.a))
    
    #---[2] transform R
    R.til<-matrix(0,nT,c(p+p1))
    for (j in 1:c(p+p1)) {
      R.j<-R[,j]
      R.j.Tn<-matrix(R.j,t,n)
      R.j.bar.a<-kronecker(colSums(R.j.Tn*w.x.a.Tn),ones.T)
      R.til[,j]<-R.j-R.j.bar.a
    }
    
    R.til.K<-t(K*R.til)
    beta.hat.x<-(chol2inv(chol(R.til.K%*%R.til)))%*%(R.til.K%*%Y.til)
    
    theta.hat[i,]<-t(beta.hat.x[1:p,])
    theta.der.hat[i,]<-t(beta.hat.x[-c(1:p),])
  }
  
  
  return(list(gradient.x=theta.hat,gradient.z=theta.der.hat,h=h,Rk=Rk))
  
})

#---[5.2] local-within local linear estimator for intercept b0(z)
np.lwll.step2<-cmpfun(function(y,x,x.eval=NULL,n,t,h,ckertype='gaussian',lono=F,Print=T){
  
  h<-as.numeric(h)
  
  y<-as.matrix(y)
  x<-as.matrix(x)
  if(is.null(x.eval)) x.eval<-x else x.eval<-as.matrix(x.eval)
  
  nT<-n*t
  nT.eval<-nrow(x.eval)
  
  d<-dim(x)[2]
  ones<-rep(1,nT)
  ones.T<-rep(1,t)
  
  
  if(ckertype=="gaussian") {
    k <- function(g) (1/sqrt(2*pi))*exp(-0.5*g^2)
  }else{
    if(ckertype=='uniform'){
      k<-function(g) ifelse(abs(g)<=1,0.5,0)
    }else{
      if(ckertype=="epanechnikov") {
        k <- function(g) ifelse(abs(g)<=1, 0.75*(1-g^2), 0)
      }else{
        if(ckertype=="biweight") {
          k <- function(g) ifelse(abs(g)<=1, (15/16)*(1-g^2)^2, 0)
        }else{
          if(ckertype=="triweight") {
            k <- function(g) ifelse(abs(g)<=1, (35/32)*(1-g^2)^3, 0)
          }
        }
      }
    }
  }
  
  der.hat<-matrix(0,nT.eval,d)
  gamma.hat<-matrix(0,nT.eval,1)
  for(i in 1:nT.eval){ 
    
    if(Print) cat(paste("\rStep 2 Estimation at",i,'of',nT.eval))
    
    K<-ones
    dx<-matrix(0,nT,d)
    for (j in 1:d){
      dx.j<-x[,j]-x.eval[i,j]
      dx[,j]<-dx.j
      v <- dx.j/h[j]
      K <- K*k(v)
    }
    
    if(lono) K[i]<-0
    
    K.Tn<-matrix(K,t,n)
    Y.Tn<-matrix(y,t,n)
    
    K.bar.n<-colSums(K.Tn)
    K.bar.nT<-kronecker(K.bar.n,ones.T)
    w.x<-K/K.bar.nT
    
    w.x.Tn<-matrix(w.x,t,n)
    y.bar<-kronecker(colSums(Y.Tn*w.x.Tn),ones.T)
    Y.til<-y-as.matrix(y.bar)
    
    X.til<-matrix(0,nT,d)
    for (j in 1:d) {
      x.j<-x[,j]
      x.j.Tn<-matrix(x.j,t,n)
      x.j.bar<-kronecker(colSums(x.j.Tn*w.x.Tn),ones.T)
      X.til[,j]<-x.j-x.j.bar
    }
    
    KX.til<-K*X.til
    X.til.K<-t(KX.til)
    beta.hat.x<-(chol2inv(chol(X.til.K%*%X.til)))%*%(X.til.K%*%Y.til)
    der.hat[i,]<-t(beta.hat.x)
    
    gamma.hat[i,]<-sum((y-dx%*%beta.hat.x)*as.matrix(w.x))/n
  }
  
  K<-ones
  for (j in 1:d){
    v <- x[,j]/h[j]
    K <- K*k(v)
  }
  
  if(lono) K[i]<-0
  
  K.Tn<-matrix(K,t,n)
  Y.Tn<-matrix(y,t,n)
  
  K.bar.n<-colSums(K.Tn)
  K.bar.nT<-kronecker(K.bar.n,ones.T)
  w.0<-K/K.bar.nT
  
  w.0.Tn<-matrix(w.0,t,n)
  y.bar<-kronecker(colSums(Y.Tn*w.0.Tn),ones.T)
  Y.til<-y-as.matrix(y.bar)
  
  X.til<-matrix(0,nT,d)
  for (j in 1:d) {
    x.j<-x[,j]
    x.j.Tn<-matrix(x.j,t,n)
    x.j.bar<-kronecker(colSums(x.j.Tn*w.0.Tn),ones.T)
    X.til[,j]<-x.j-x.j.bar
  }
  
  KX.til<-K*X.til
  X.til.K<-t(KX.til)
  beta.hat.0<-(chol2inv(chol(X.til.K%*%X.til)))%*%(X.til.K%*%Y.til)
  gamma.hat.0<-sum((y-x%*%beta.hat.0)*as.matrix(w.0))/n
  
  m.hat<-gamma.hat-gamma.hat.0
  
  return(list(m.hat=m.hat,gradient.x=der.hat,h=h))
  
})

#---[6.1] local-within local linear estimator for one-way VCM: X'b(Z)
#   Note: it assumes b0(z)=0, so only returns function and derivative of b1(z),...,bp(z)
vcm.nc.lwll<-cmpfun(function(y,x,z,z.eval=NULL,n,t,C.h=NULL,ckertype='gaussian'){
  
  y= as.matrix(y); x= as.matrix(x);  z= as.matrix(z);  nT=n*t
  
  if(is.null(z.eval)) z.eval=z else z.eval=as.matrix(z.eval)
  nT.eval<-nrow(z.eval)
  p<- ncol(x)
  q<- ncol(z)
  
  ones<-rep(1,nT)
  ones.T<-rep(1,t)
  
  if(ckertype=="gaussian") {
    k <- function(g) (1/sqrt(2*pi))*exp(-0.5*g^2)
    if(is.null(C.h)==TRUE){
      h<-1.059*apply(z,2,sd)*nT^(-1/(4+q))
    }else{
      h<-C.h*apply(z,2,sd)*nT^(-1/(4+q))
    }
  }else{
    if(ckertype=='uniform'){
      k<-function(g) ifelse(abs(g)<=1,0.5,0)
      if(is.null(C.h)==TRUE){
        h<-1.843*apply(z,2,sd)*nT^(-1/(4+q))
      }else{
        h<-C.h*apply(z,2,sd)*nT^(-1/(4+q))
      }
    }else{
      if(ckertype=="epanechnikov") {
        k <- function(g) ifelse(abs(g)<=1, 0.75*(1-g^2), 0)
        if(is.null(C.h)==TRUE){
          h<-2.345*apply(z,2,sd)*nT^(-1/(4+q))
        }else{
          h<-C.h*apply(z,2,sd)*nT^(-1/(4+q))
        }
      }else{
        if(ckertype=="biweight") {
          k <- function(g) ifelse(abs(g)<=1, (15/16)*(1-g^2)^2, 0)
          if(is.null(C.h)==TRUE){
            h<-2.778*apply(z,2,sd)*nT^(-1/(4+q))
          }else{
            h<-C.h*apply(z,2,sd)*nT^(-1/(4+q))
          }
        }else{
          if(ckertype=="triweight") {
            k <- function(g) ifelse(abs(g)<=1, (35/32)*(1-g^2)^3, 0)
            if(is.null(C.h)==TRUE){
              h<-3.154*apply(z,2,sd)*nT^(-1/(4+q))
            }else{
              h<-C.h*apply(z,2,sd)*nT^(-1/(4+q))
            }
          }
        }
      }
    }
  }
  theta.hat<-matrix(0,nT.eval,p)
  theta.der.hat<-matrix(0,nT.eval,p*q)
  for(i in 1:nT.eval){ 
    
    R<-x
    K<-ones
    for (jj in 1:q){
      dz<-(z[,jj]-z.eval[i,jj])
      v.z<-dz/h[jj]
      K<-K*k(v.z)
      
      for (jjj in 1:p){
        x.dz <- x[,jjj]*dz
        R <- cbind(R,x.dz)
      }
    }
    
    #---[a]
    K.Tn<-matrix(K,t,n)
    
    K.bar.n<-colSums(K.Tn)
    K.bar.nT<-kronecker(K.bar.n,ones.T)
    w.x.a<-K/K.bar.nT
    w.x.a.Tn<-matrix(w.x.a,t,n)
    
    #---[1] transform Y
    Y.Tn<-matrix(y,t,n)
    y.bar.a<-kronecker(colSums(Y.Tn*w.x.a.Tn),ones.T)
    
    Y.til<-y-as.matrix((y.bar.a))
    
    #---[2] transform R
    R.til<-matrix(0,nT,2*p)
    for (j in 1:c(2*p)) {
      R.j<-R[,j]
      R.j.Tn<-matrix(R.j,t,n)
      R.j.bar.a<-kronecker(colSums(R.j.Tn*w.x.a.Tn),ones.T)
      R.til[,j]<-R.j-R.j.bar.a
    }
    
    R.til.K<-t(K*R.til)
    beta.hat.x<-(chol2inv(chol(R.til.K%*%R.til)))%*%(R.til.K%*%Y.til)
    
    theta.hat[i,]<-t(beta.hat.x[1:p,])
    theta.der.hat[i,]<-t(beta.hat.x[-c(1:p),])
  }
  
  return(list(gradient.x=theta.hat,gradient.z=theta.der.hat,h=h))
  
})

#---[6.2] Objective function for CVLS bandwidth
lwll.cvls<-function(h,y,x,z,ckertype){
  
  h<-as.numeric(h)
  
  vcm.nc<-vcm.lwll.step1(y=y,x=x,z=z,z.eval=z,n=n,t=t,h=h,ckertype = ckertype,lono=T,Print = F)
  b.hat<-vcm.nc$gradient.x
  y.hat<-y-rowSums(x*b.hat)
  
  npm<-np.lwll.step2(y=y.hat,x=z,x.eval=z,n=n,t=t,h=h,ckertype = ckertype,lono=T,Print = F)
  m.hat<-npm$m.hat
  
  e.hat<-y.hat-m.hat
  e.til.hat<-within.tran(e.hat,n,t)
  
  return(mean((e.til.hat)^2))
  
}

#==============================================================================================================================
#---[7] YWC local-within local linear VCM: b0(z)+x'b(z): using [1]-[6.2]
#---y: dependent variable
#---x: regressors with smooth coefficients (of dimension p)
#---z: smooth variables (of dimention q)
#---z.eval: evaluation points at which functions are estimated
#---n,t: sample size of individuals and time, respectively
#---bwmethod: bandwidth selection method, including rule-of-thumb ('rot'), cross-validation ('cv'), and given by users ('given')
#---C.h: the scaling factor of bandwidth, which is 1.059 by default.
#---eps.hat: residual
#---eval: an indicator for whether evaluation data is used. If eval=TRUE (by default), then z.eval is not equal to z
#---eval.type: the type of evaluation points, including 'sparse' (a subsample of z), 'grid' (a grid point from z's support), or
#   'given', allowing for a particular choice by the users
#---grid.length: the total evaluation points if eval.type='grid
#---sparse.by is the distance between each of the chosen sample points. Use this when specifying 
#   eval.type='sparse'. Default is 1
#---ckertype: contineous kernel type. The default is gaussian kernel

YWC.vcm<-cmpfun(function(y,x,z,z.eval=NULL,n,t,bwmethod=c('rot','cv','given'),C.h=1.059,
                              eps.hat=NULL,eval=TRUE,eval.type=c('sparse','grid','given'),grid.length=50,sparse.by=1,ckertype='gaussian'){
  
  #---Data checking
  y= as.matrix(y); x= as.matrix(x);  z= as.matrix(z)
  ny<-nrow(y); nz<-nrow(z); nx<-nrow(x)
  
  if(ny!=nz | ny!=nx | nz!=nx) stop("At least one of variables from (y,x,z) contains different observations")
  
  nT=n*t; p<- ncol(x); q<- ncol(z); p1<-p+1
  x.til<-cbind(rep(1,nT),x)
  
  #---Evaluation points
  if(eval){
    if(eval.type=='sparse'){
      if(sparse.by>=nT) stop('Evaluation point is only one. Reset sparse.by<=nT')
      if(sparse.by<1)  stop('Evaluation points are higher than original points and thus too slow. Reset sparse.by>1')
      
      eval.index<-seq(1,nT,by=sparse.by)
      z.eval<-z[eval.index,,drop=FALSE]
      nT.eval<-nrow(z.eval)
    }else{
      if(eval.type=='grid'){
        z.eval<-sapply(1:q, function(s){seq(quantile(z[,s],0.01),quantile(z[,s],0.99),length=grid.length)})
        nT.eval<-nrow(z.eval)
      }else{
        if(eval.type=='given'){
          if(is.null(z.eval)) stop("eval.type='given' requires your specification of evaluation points")
          z.eval<-as.matrix(z.eval);
          if(ncol(z.eval)!=q) stop("Provided z.eval matrix does not have the same column size to z")
          nT.eval<-nrow(z.eval)
        }
      }
    }
  }else{
    z.eval<-z
  }
  
  nT.eval<-nrow(z.eval)
  
  #---Bandwidth selection
  
  if(bwmethod=='rot'){
    cat(paste("Bandwidth selected based on Rule of Thumb with scaling factor",C.h,'\n'))
    h<-C.h*apply(z,2,sd)*nT^(-1/(4+q))
    cat("\n")
  }else{
    if(bwmethod=='cv'){
      cat(paste("Bandwidth selected based on CVLS. Computing...",'\n'))
      cat("\n")
      bw.start <- C.h*apply(z,2,sd)*(nT^(-1/(4+q)))
      lower<-0; upper<-10*sd(z)
      
      bw.optim <- bobyqa(bw.start,lwll.cvls,lower,upper,y=y,x=x,z=z,ckertype=ckertype)
      h<-bw.optim$par
      error.code<-bw.optim$ierr
      
      cat(paste("Optimal CV Bandwidth is found by",round(h,4),"with error code",error.code,'\n'))
      cat("\n")
    }else{
      if(bwmethod=='given'){
        cat(paste("Bandwidth selected based on author's choice",'\n'))
        cat("\n")
        h<-h 
        if(length(which(h<=0))!=0) stop("Provided h is not positive")
      }
    }
  }
  
  #---Perform estimation
  
  
  cat(paste("Step 1 Estimation...",'\n'))
  vcm.nc<-vcm.lwll.step1(y=y,x=x,z=z,z.eval=z.eval,n=n,t=t,h=h,ckertype = ckertype)
  bx.hat<-vcm.nc$gradient.x
  b.der.hat<-vcm.nc$gradient.z
  h.step1<-vcm.nc$h
  Rk<-vcm.nc$Rk
  
  cat('\n')
  cat(paste("Step 2 Estimation...",'\n'))
  vcm.nc.full<-vcm.lwll.step1(y=y,x=x,z=z,z.eval=z,n=n,t=t,h=h,ckertype = ckertype)
  b.hat.full<-vcm.nc.full$gradient.x
  xb.hat<-rowSums(x*b.hat.full)
  y.hat<-y-xb.hat
  
  npm<-np.lwll.step2(y=y.hat,x=z,x.eval=z.eval,n=n,t=t,h=h,ckertype = ckertype)
  m.hat<-npm$m.hat
  h.step2<-npm$h
  
  b.hat<-cbind(m.hat,bx.hat)
  
  
  if(is.null(eps.hat)){
    
    cat('\n')
    cat('\n')
    cat(paste("Constructing inference of step 1 and 2 estimates...",'\n'))
    
    if(nT!=nT.eval){
      
      cat(paste("Re-performing step 2 using all points of z to construct residual",'\n'))
      npm.full<-np.lwll.step2(y=y.hat,x=z,x.eval=z,n=n,t=t,h=h,ckertype = ckertype)
      e.hat<-y.hat-npm.full$m.hat
      eps.hat<-as.numeric(within.tran(e.hat,n,t,constrain = F))
    }else{
      e.hat<-y.hat-m.hat
      eps.hat<-as.numeric(within.tran(e.hat,n,t,constrain = F))
    }
  }else{
    cat(paste("Eps.hat is given by the author",'\n'))
    eps.hat<-as.numeric(eps.hat)
    if(nrow(eps.hat)!=nT) stop("The provided eps.hat is different from nT")
  }
  
  fz.hat<-npuden(x=z,x.eval=z.eval,ckertype = ckertype)$fhat
  eps.hat.sq<-eps.hat^2
  sigma.eq.z<-lcls(y=eps.hat.sq,x=z,x.eval = z.eval,ckertype = ckertype)$fitted.y
  
  #---XX' for upper-left block matrix
  XXi.box<-array(0,dim=c(p,p,nT))
  for (ii in 1:p) {
    XX.tr<-as.matrix(x[,ii]*x[,c(ii:p)])
    for (jj in ii:p) {
      XXi.box[ii,jj,]=XXi.box[jj,ii,]<-XX.tr[,c(jj-(ii-1))]
    }
  }
  
  #---X.tilX.til' for lower-right block matrix
  XXi.tilde.box<-array(0,dim=c(p1,p1,nT))
  for (ii in 1:p1) {
    XX.til.tr<-as.matrix(x.til[,ii]*x.til[,c(ii:p1)])
    for (jj in ii:p1) {
      XXi.tilde.box[ii,jj,]=XXi.tilde.box[jj,ii,]<-XX.til.tr[,c(jj-(ii-1))]
    }
  }
  
  #---[1.1] step 1 estimator inference: S
  
  EX.EXtr.box<-array(0,dim=c(p,p,nT.eval))
  EX.Z<-sapply(1:p, function(s){llls(y=x[,s],x=z,x.eval=z.eval,ckertype = ckertype)$fitted.y})
  for (ii in 1:p) {
    EX.ZEX.Z.tr<-as.matrix(EX.Z[,ii]*EX.Z[,c(ii:p)])
    for (jj in ii:p) {
      EX.EXtr.box[ii,jj,]=EX.EXtr.box[jj,ii,]<-EX.ZEX.Z.tr[,c(jj-(ii-1))]
    }
  }
  
  
  #---upper-left block matrix
  Omega.bar<-array(0,dim=c(p,p,nT.eval))
  for (ii in 1:p) {
    
    EXX.tr.grid<-sapply(ii:p, function(s){llls(y=XXi.box[ii,s,],x=z,x.eval=z.eval,ckertype = ckertype,C.h = C.h)$fitted.y})
    EXEX.tr.grid<-sapply(ii:p, function(s){EX.EXtr.box[ii,s,]})
    
    S1.z<-EXX.tr.grid-EXEX.tr.grid
    
    for (jj in ii:p) {
      Omega.bar[ii,jj,]=Omega.bar[jj,ii,]<-S1.z[,c(jj-(ii-1))]
    }
  }
  
  
  #---lower-right block matrix
  omega<-array(1,dim=c(p1,p1,nT.eval))
  mu.2.k<-diag(1,q,q)
  for (ii in 1:p1) {
    
    S2.z<-sapply(ii:p1, function(s){llls(y=XXi.tilde.box[ii,s,],x=z,x.eval=z.eval,ckertype = ckertype,C.h = C.h)$fitted.y})
    
    for (jj in ii:p1) {
      omega[ii,jj,]=omega[jj,ii,]<-S2.z[,c(jj-(ii-1))]
    }
  }
  
  Omega<-array(1,dim=c(p1*q,p1*q,nT.eval))
  for (i in 1:nT.eval) {
    Omega[,,i]<-kronecker(omega[,,i],mu.2.k)
  }
  
  S<-array(0,dim=c(c(p+p1*q),c(p+p1*q),nT.eval))
  S[c(1:p),c(1:p),]<-Omega.bar
  S[c(p1:c(p+p1*q)),c(p1:c(p+p1*q)),]<-Omega
  
  
  #---[1.2] step 1 estimator inference: Omega
  
  #---upper-left block matrix
  Sigma.F<-array(0,dim=c(p,p,nT.eval))
  for (i in 1:nT.eval) {
    Sigma.F[,,i]<-(Rk^q)*sigma.eq.z[i]*Omega.bar[,,i]
  }
  
  #---lower-right block matrix
  R2.k<-diag(1/(4*sqrt(pi)),q,q)
  Sigma.D<-array(1,dim=c(p1*q,p1*q,nT.eval))
  for (i in 1:nT.eval) {
    Sigma.D[,,i]<-kronecker(omega[,,i]*sigma.eq.z[i],R2.k)
  }
  
  Sigma<-array(0,dim=c(c(p+p1*q),c(p+p1*q),nT.eval))
  Sigma[c(1:p),c(1:p),]<-Sigma.F
  Sigma[c(p1:c(p+p1*q)),c(p1:c(p+p1*q)),]<-Sigma.D
  
  Dh.inv<-diag(0,p+p1*q,p+p1*q)
  diag(Dh.inv)<-c(rep(1,p),rep(1/h,p1*q))
  h.prod<-prod(h)
  
  Asy.theta_1.sd<-matrix(0,nT.eval,c(p+p1*q))
  for (iii in 1:nT.eval){
    S.inv.i<-solve(S[,,iii])
    Omega.i<-Sigma[,,iii]
    S.inv_O_S.inv_fhat<-((1/(nT*h.prod*fz.hat[iii]))*Dh.inv)%*%S.inv.i%*%Omega.i%*%S.inv.i
    Asy.theta_1.sd[iii,]<-sqrt(abs(diag(S.inv_O_S.inv_fhat)))
  }
  
  #---[2] step 2 estimator inference
  
  #---construct Vz
  V.z<-matrix(0,nT.eval,1)
  for (i in 1:nT.eval) {
    
    Omega.bar.inv<-solve(Omega.bar[,,i])
    
    x_Ex.z_i<-sapply(1:p,function(s){x[,s]-EX.Z[i,s]})
    v.z.inner.i<-t(EX.Z[i,,drop=F]%*%Omega.bar.inv%*%t(x_Ex.z_i))
    Ev.z.inner.i<-llls(y=v.z.inner.i^2,x=z,x.eval=z.eval[i],ckertype = ckertype,C.h = C.h)$fitted.y
    
    V.z[i,]<-(1/fz.hat[i])*sigma.eq.z[i]*(1/sqrt(8*pi))*( (1/sqrt(4*pi))^q + Ev.z.inner.i )
    
  }
  
  #---construct V0
  f0.hat<-npuden(x=z,x.eval=0,ckertype = ckertype)$fhat
  sigma.eq.0<-lcls(y=eps.hat.sq,x=z,x.eval = 0,ckertype = ckertype)$fitted.y
  
  EX.EXtr.0<-array(0,dim=c(p,p,1))
  EX.0<-sapply(1:p, function(s){llls(y=x[,s],x=z,x.eval=0,ckertype = ckertype)$fitted.y})
  for (ii in 1:p) {
    EX.ZEX.Z.tr<-EX.0[ii]*EX.0[c(ii:p)]
    for (jj in ii:p) {
      EX.EXtr.0[ii,jj,]=EX.EXtr.0[jj,ii,]<-EX.ZEX.Z.tr[c(jj-(ii-1))]
    }
  }
  
  Omega.bar.0<-array(0,dim=c(p,p,1))
  for (ii in 1:p) {
    
    EXX.tr.0<-sapply(ii:p, function(s){llls(y=XXi.box[ii,s,],x=z,x.eval=0,ckertype = ckertype,C.h = C.h)$fitted.y})
    EXEX.tr.0<-sapply(ii:p, function(s){EX.EXtr.0[ii,s,]})
    
    S1.0<-EXX.tr.0-EXEX.tr.0
    
    for (jj in ii:p) {
      Omega.bar.0[ii,jj,]=Omega.bar.0[jj,ii,]<-S1.0[c(jj-(ii-1))]
    }
  }
  
  
  Omega.bar.0.inv<-solve(Omega.bar.0[,,1])
  
  x_Ex.0<-sapply(1:p,function(s){x[,s]-EX.0[s]})
  v.z.inner.0<-t(t(as.matrix(EX.0))%*%Omega.bar.inv%*%t(x_Ex.0))
  Ev.z.inner.0<-llls(y=v.z.inner.0^2,x=z,x.eval=0,ckertype = ckertype,C.h = C.h)$fitted.y
  
  V.0<-(1/f0.hat)*sigma.eq.0*(1/sqrt(8*pi))*( (1/sqrt(4*pi))^q + Ev.z.inner.0 )
  
  Asy.b0.sd<-sqrt((1/(nT*h.prod))*(V.z+V.0))
  
  beta.hat.sd<-cbind(Asy.b0.sd,Asy.theta_1.sd[,c(1:p)])
  beta.der.hat.sd<-cbind(Asy.theta_1.sd[,-c(1:p)])
  
  cat("\n")
  cat("\n")
  cat(paste("Model Estimation Complete",'\n'))
  
  return(list(y=y,x=x,z=z,z.eval=z.eval,n=n,t=t,eval.type=eval.type,grid.length=grid.length,h=h,
              beta.hat=b.hat,beta.hat.sd=beta.hat.sd,beta.der.hat=b.der.hat,beta.der.hat.sd=beta.der.hat.sd,
              h.step1=h.step1,h.step2=h.step2,bwmethod=bwmethod,ckertype=ckertype))
})
#==============================================================================================================================

#==============================================================================================================================
#---[8] Results for YWC local-within local linear VCM
#---ywc.model: the estimation from YWC.vcm.lwll
#---alpha.table is a set of chosen percentiles of the distribution of function estimates. 
#   For example, alpha.table=c(0.25, 0.5, 0.75) means that the numerical estimation results of frontier
#   functions will be reported at its 25-th, 50-th, and 75-th percentile associated with its standard error. 
#---figure is a logical value to indicate whether the estimated functions are ploted. Press <ENTER> in 
#   console to see each plot. Default is TRUE.
#---margin.b is the a by b size of figure frame for plotting function estimates \phi(z) and b(z)

result.YWC.vcm<-function(ywc.model,alpha.table=NULL,figure=TRUE,alpha.CI=NULL,margin.b=c(2,2)){
  
  y=ywc.model$ywc.model$y;x=ywc.model$x;z=ywc.model$z;z.eval=ywc.model$z.eval;
  p<-ncol(x); q<-ncol(z); p1<-p+1
  
  n=ywc.model$n;t=ywc.model$t; nT=n*t; nT.eval<-nrow(z.eval); h<-ywc.model$h
  
  eval.type=ywc.model$eval.type;grid.length=ywc.model$grid.length;
  
  beta.hat=ywc.model$beta.hat; beta.hat.sd=ywc.model$beta.hat.sd;
  beta.der.hat=ywc.model$beta.der.hat; beta.der.hat.sd=ywc.model$beta.der.hat.sd;
  
  h.step1=ywc.model$h.step1;h.step2=ywc.model$h.step2;bwmethod=ywc.model$bwmethod; ckertype=ywc.model$ckertype
  
  y.name<-if(is.null(colnames(y))) "Y" else colnames(y)
  z.name<-if(is.null(colnames(z))) paste0("Z",c(1:q)) else colnames(z)
  x.name<-if(is.null(colnames(x))) paste0("X",c(1:p)) else colnames(x)
  inter.fun.name<-paste('b0(',z.name,')')
  
  
  #---[1] summary of individual function and derivative
  if(is.null(alpha.table)){
    Q=c(0.1, 0.25, 0.5, 0.75, 0.9); LQ=length(Q)
  }else{
    if(alpha.table>1 | alpha.table<0) stop("'alpha.table' must be chosen between 0 and 1 for quantile of variables in Tables " )
    Q<-seq(0,1,alpha.table)
    Q<-Q[-c(1,length(Q))]; Q=c(0.1, Q, 0.9)
    LQ<-length(Q)
  }
  
  row.name<-rep(NA,2*LQ)
  Q.ch<-paste0(Q)
  for (s in 1:LQ){row.name[c((2*s-1):(2*s))]<-c(Q.ch[s],'')}
  
  M.m<-matrix(0,LQ*2,p1)
  M.m.der<-matrix(0,LQ*2,p1)
  for (d in 1:p1) {
    
    G.m<-cbind(z.eval,beta.hat[,d],beta.hat.sd[,d])
    G.m.der<-cbind(z.eval,beta.der.hat[,d],beta.der.hat.sd[,d])
    
    for (jjj in 1:LQ){
      cond.fun.m<-t(G.m[which.min(abs(G.m[,1]-quantile(G.m[,1],Q[jjj]))),c(-1)])
      M.m[c(2*jjj-1):(2*jjj),d]<-round(cond.fun.m,4)
      
      cond.fun.der<-t(G.m.der[which.min(abs(G.m.der[,1]-quantile(G.m.der[,1],Q[jjj]))),c(-1)])
      M.m.der[c(2*jjj-1):(2*jjj),d]<-round(cond.fun.der,4)
    }
  }
  
  M.m<-round(rbind(colMeans(beta.hat),M.m),4)
  colnames(M.m)<-paste0('b',c(1:p1),'(',z.name,')')
  rownames(M.m)<-c('mean',row.name)
  
  M.m.der<-round(rbind(colMeans(beta.der.hat),M.m.der),4)
  colnames(M.m.der)<-paste0("b",c(1:p1),"'(",z.name,')')
  rownames(M.m.der)<-c('mean',row.name)
  
  cat("\n")
  cat('============================================================================================','\n')
  cat('\n')
  cat("A FE Local-Within Varying Coefficient Model by Yao, Wang, Cai (2024)             \n")
  cat('\n')
  cat('Dependent variables:                            '     ,y.name,'\n')
  cat('Smooth variables:                               '     ,z.name,'\n')
  cat('Linear variables:                               '     ,x.name,'\n')
  cat('\n')
  cat('Summary of Smooth Coefficient function bj(.) conditional on quantile of Z'          ,'\n')
  cat('--------------------------------------------------------------------------------------------','\n')
  cat('\n')
  cat('[A.1]: Function bj(.) Summary:'                                                              ,'\n')
  cat('\n')
  print(M.m)
  cat('\n')
  cat('[A.2]: Derivative bj(.) Summary:'                                                            ,'\n')
  cat('\n')
  print(M.m.der)
  cat('\n')
  cat('--------------------------------------------------------------------------------------------')
  cat('\n')
  cat("Firms:                         ", n,           "\n")
  cat("Time Span:                     ", t,           "\n")
  cat("Total Obs:                     ", nT,          "\n")
  cat("Total Eval:                    ", nT.eval,     "\n")
  cat("Bandwith method:               ", bwmethod,    "\n")
  cat("Bandwidth used:                ", round(h,4),  "\n")
  cat('============================================================================================','\n')
  cat('\n')
  
  if(figure){
    
    if(is.null(alpha.CI)){
      alpha.CI=-qnorm((0.05/2)); a=0.05
    }else{
      alpha.CI=-qnorm((alpha.CI/2)); a=alpha.CI
    }
    
    par(mfrow=margin.b)  
    
    z.range<-range(z.eval)
    b.beta.CI<-array(0,dim=c(nT.eval,4,p1))
    for (d in 1:p1) {
      GM<-cbind(z.eval,beta.hat[,d],beta.hat[,d]-alpha.CI*beta.hat.sd[,d],beta.hat[,d]+alpha.CI*beta.hat.sd[,d])
      b.beta.CI[,,d]<-GM[order(GM[,1]),]
    }
    
    for (sss in 1:p1) {
      y.range<-c(min(b.beta.CI[,3,sss]),max(b.beta.CI[,4,sss]))
      plot(b.beta.CI[,1,sss],b.beta.CI[,2,sss],type='l',col='red',xlab=z.name,xlim=z.range
           ,ylab=paste0("Estimated function of b",sss,'(',z.name[sss],')'),lwd=2,ylim=y.range
           ,main=paste0("b",sss,'(',z.name,')'))
      lines(b.beta.CI[,1,sss],b.beta.CI[,3,sss],type='l',col='black',lty=2,lwd=2)
      lines(b.beta.CI[,1,sss],b.beta.CI[,4,sss],type='l',col='black',lty=2,lwd=2)
      legend('topleft',c(paste0('b',sss,'(',z.name,')'),paste0(round((1-a)*100)," %CI")),lwd=c(2),lty=c(1,2),col=c('red','black'))
      abline(h=0, lty=3) 
    }
    
    readline(cat ("Press <Enter> to continue for derivative estimates of b'(Z): "))
    par(mfrow=margin.b)
    b.beta.der.CI<-array(0,dim=c(nT.eval,4,p1))
    for (d in 1:p1) {
      GM<-cbind(z.eval,beta.der.hat[,d],beta.der.hat[,d]-alpha.CI*beta.der.hat.sd[,d],beta.der.hat[,d]+alpha.CI*beta.der.hat.sd[,d])
      b.beta.der.CI[,,d]<-GM[order(GM[,1]),]
    }
    
    for (sss in 1:p1) {
      y.range<-c(min(b.beta.der.CI[,3,sss]),max(b.beta.der.CI[,4,sss]))
      plot(b.beta.der.CI[,1,sss],b.beta.der.CI[,2,sss],type='l',col='red',xlab=z.name,xlim=z.range
           ,ylab=paste0("Estimated derivative of b'",sss,'(',z.name[sss],')'),lwd=2,ylim=y.range
           ,main=paste0("b'",sss,'(',z.name,')'))
      lines(b.beta.der.CI[,1,sss],b.beta.der.CI[,3,sss],type='l',col='black',lty=2,lwd=2)
      lines(b.beta.der.CI[,1,sss],b.beta.der.CI[,4,sss],type='l',col='black',lty=2,lwd=2)
      legend('topleft',c(paste0('b',sss,'(',z.name,')'),paste0(round((1-a)*100)," %CI")),lwd=c(2),lty=c(1,2),col=c('red','black'))
      abline(h=0, lty=3) 
    }
    
  }
  result.list<-list(b.beta.CI=b.beta.CI,b.beta.der.CI=b.beta.der.CI)
}
#==============================================================================================================================



#================================
#---[Ex] q=1
#================================

n=100;t=10; p=1

Dgp<-DGP(dgp.type = 2,p=p, grid.length=50,n=n,t=t,fe.normalize = F,a=-2,b=2)

y= as.matrix(Dgp$y); x= as.matrix(Dgp$x);  z= as.matrix(Dgp$z); z.eval=as.matrix(Dgp$z.eval)

ywc.model<-YWC.vcm(y=y,x=x,z=z,z.eval=z.eval,n=n,t=t,bwmethod='cv',eval.type = 'given')
result.model<-result.YWC.vcm(ywc.model,margin.b=c(1,2))