Quantnet :: Show

rm(list=ls(all=TRUE))
library(KernSmooth)
library(VGAM)
library(expectreg)
library(mboost)
library(KernSmooth) #install.packages("KernSmooth") if package missing
lprq2 <- function(x, y, h, tau, x0) # modified from lprq, s.t. we can specify where to estimate quantiles
        {xx <- x0
        fv <- xx
        dv <- xx
        for(i in 1:length(xx)) {
                z <- x - xx[i]
                wx <- dnorm(z/h)
                r <- rq(y~z, weights=wx, tau=tau, ci=FALSE)
                fv[i] <- r$coef[1]
                dv[i] <- r$coef[2]
                               }    
        list(xx = xx, fv = fv, dv = dv)
        }
lprq3 <- function(x, y, h, x0) # modified from lprq, s.t. we can specify where to estimate quantiles, but random quantiles
       {xx <- x0
        fv <- xx
        dv <- xx
        for(i in 1:length(xx)) {
                z <- x - xx[i]
                wx <- dnorm(z/h)
                r <- rq(y~z, weights=wx, tau=runif(1), ci=FALSE)
                fv[i] <- r$coef[1]
                dv[i] <- r$coef[2]
                               }  
        list(xx = xx, fv = fv, dv = dv)
        }
local.expectile=function(y,x,z,h,q){
   # parameters: x=variable; h=bandwidth; z=grid point; ker=kernel
   nz<-length(z)
   nx<-length(x)
   x_0=rep(1,nx*nz)
   dim(x_0)=c(nx,nz)
   x_1=t(x_0)
   x0=x*x_0 
   w1 <- 0 * y + 0.5 
   y0=y*x_0 
   #y1=0.5*y
   y11<-0.5*y*x_1
   x1=z*x_1
   x10=x0-t(x1)
    #if(ker==1){x1=kernel(x0/h)}         # Epanechnikov kernel
    x11=kernel(x10/h)   
   #x2=(y1<y0)*x_1
   #for (i in 1:nz)
 it=1 
 dw1 = 1
 while (dw1!=0 & it < 100)
  #while ( it < 500)
 { 
   v1=((y0<t(y11))-2*q*(y0<t(y11))+q)*x11
   v2=y*v1
   f1=apply(v1,2,mean)
   f2=apply(v2,2,mean)
   f3=f2/f1
   w01 <- w1
   w1 <- as.vector(ifelse(y0>t(f3*x_1), q, 1 - q))    
   dw1 <- sum(w1 != w01, na.rm = TRUE)
   y1=f3
   y11=f3*x_1
   it=it+1
   }
   return(y1)
}
kernelq <- function(u){dnorm(u,mean=0,sd=1)} 
kernel<-function(x){0.75*(1-x^2)*(abs(x)<=1)}
coverate<-0
area.cover<-0
B<-500
n=500
q=0.95
h=0.2
alpha=0.05
gridn<-n
cc    <- 1/4
set.seed(10*3)
x<-runif(n,0,2)
y<-1.5*x+2*sin(pi*x)+rnorm(n)
e_theor=1.5*x+2*sin(pi*x)+enorm(q)
z<-x
v   <- x # v for the nonparametric part
bound<- c(0, 1)
yuv  <-sort(v) # just sort them for later use
yur  <-y[order(v)]
yuv  <-sort(v) # just sort them for later use
yur  <-y[order(v)]
h12   <-(0.5*(1-0.5)/dnorm(qnorm(0.5))^2)^0.2*h
b2    <-max(h/sd(yuv)*sd(yur), h12^5/(h)^3, h/10)
lambda<- 1/(2*sqrt(pi)) # this is for normal kernel, if quartic kernel, value is 5/7
delta<--log(h)/log(n)
dd    <- sqrt(2*delta*log(n))+(2*delta*log(n))^(-1/2)*log(cc/2/pi)
fxd   <-bkde(yuv, gridsize = gridn, range.x = c(min(z), max(z)))
values<-local.expectile(y,x,z,h,q)
yy1<-y-values[order(x)]
fx_yy  <-bkde(sort(yy1), gridsize = gridn, range.x = c(min(yy1), 0))
delta_yy<-c(min(yy1),fx_yy$x)
sigma<-q^2*mean(yy1^2)+(1-2*q)*sum(yy1^2*fx_yy$y*diff(delta_yy))
f_x<-ecdf(yy1)
fx<-q+(1-2*q)*f_x(0)
bandt <-(fxd$y)^(1/2)*fx
cn    <-log(2)-log(abs(log(1-alpha)))
band  <-(n*h)^(-1/2)*sqrt(lambda*sigma)*bandt^(-1)*(dd+cn*(2*delta*log(n))^(-1/2))
coverate1<-min(e_theor[order(x)]<=values[order(x)]+band,e_theor[order(x)]>=values[order(x)]-band)
coverate<-coverate+coverate1
area.cover<-mean(band)+area.cover
plot(x,y, type="p", pch=20,  ylim=c(min(y),max(y)),xlab="X", ylab="Y", cex.lab=1.8, cex.axis=1.8)
lines(sort(x),values[order(x)],col="blue", lwd=2, type="l")
lines(sort(x),e_theor[order(x)], lwd=2, col=1)
lines(sort(x), values[order(x)]-band,col = "red", lty = 2, lwd=3)
lines(sort(x), values[order(x)]+band,col = "red", lty = 2, lwd=3)
win.graph()
library(VGAM)
library(expectreg)
x<-seq(0.0001,0.9999,length=1000)
y1<-qenorm(x)
y2<-qnorm(x)
y<-cbind(y1,y2,0)
matplot(x,y,type="l", lty=1, cex.axis=2, cex.lab = 2, xlab="tau", ylab= "", lwd=3,col=c(3,4,2),ylim=c(-2,2))