solution week 2.R

richard — Sep 24, 2013, 2:39 PM

c <- 5
p <- 0.1
n <- 1000
z <- rnorm(n)
e <- ifelse((runif(n)>p),z, c*z)
hist(e,freq=F,breaks="FD",ylim=c(0,0.4),xlim=c(-15,15))
xp <- seq(-15,15,length=1000)
lines(xp,(1-p)*dnorm(xp)+p*dnorm(xp,0,c),col="blue")
lines(density(e),col="red")

# p<-0.1
# c <- 1.5

beta <- 1
n <- 20
x <- runif(n)
z <- rnorm(n)
e <- ifelse(runif(n)>p,z, c*z)
y <- beta*x+e

plot(x,y)
abline(c(0,beta))
abline(lm(y~x-1),col="red")

plot(x,y,xlim=c(-4,5),ylim=c(-4,5))
abline(c(0,beta))
abline(lm(y~x-1),col="red")

c <- 1.5
p <- 0.1

beta <- 1
n <- 20
c <- 5
p <- 0.1

betahatLS <- numeric(1000)
sigmahatLS <- numeric(1000)
for (i in 1:1000) {
    x <- runif(n)
    z <- rnorm(n)
    e <- ifelse(runif(n)>p,z, c*z)
    y <- beta*x+e
    fitLS <- lm(y~x-1)
    betahatLS[i] <- fitLS$coef
    sigmahatLS[i] <- summary(fitLS)$sigma
    }

library(quantreg)

Loading required package: SparseM

Attaching package: 'SparseM'

The following object is masked from 'package:base':

backsolve

betahatRQ <- numeric(1000)
sigmahatRQ <- numeric(1000)
for (i in 1:1000) {
    x <- runif(n)
    z <- rnorm(n)
    e <- ifelse(runif(n)>p,z, c*z)
    y <- beta*x+e
    fitRQ <- rq(y~x-1)
    betahatRQ[i] <- fitRQ$coef
    sigmahatRQ[i] <- sqrt( sum((fitRQ$res)^2) / (n-1) )
    }


p <- 0
beta <- 1
n <- 20
c <- 2

betahatRQ <- numeric(1000)
sigmahatRQ <- numeric(1000)
sigmahatRQalt <- numeric(1000)
for (i in 1:1000) {
    x <- runif(n)
    z <- rnorm(n)
    e <- ifelse(runif(n)>p,z, c*z)
    y <- beta*x+e
    fitRQ <- rq(y~x-1)
    betahatRQ[i] <- fitRQ$coef
    sigmahatRQ[i] <- sqrt( sum((fitRQ$res)^2) / (n-1) )
    sigmahatRQalt[i] <-  sum(abs(fitRQ$res)) / n
    }

sd(betahatRQ)

[1] 0.5092

mean(betahatRQ)-beta

[1] 0.0007067

sqrt(mean((betahatRQ-beta)^2))

[1] 0.5089

sigma <- 1


sd(sigmahatRQ)

[1] 0.171

mean(sigmahatRQ)-sigma

[1] 0.0003926

sqrt(mean((sigmahatRQ-sigma)^2))

[1] 0.1709

sd(sigmahatRQalt)

[1] 0.1363

mean(sigmahatRQalt)-sigma

[1] -0.2356

sqrt(mean((sigmahatRQalt-sigma)^2))

[1] 0.2722

mean(abs(rnorm(1000)))

[1] 0.7851

sqrt(mean(rnorm(1000)^2))

[1] 0.9781

sqrt(2/pi)

[1] 0.7979

N <- 1000000
mean(abs(rnorm(N)))

[1] 0.7982

mean(rnorm(N)^2)

[1] 1.001

p <- 0
c <- 2
beta <- 0
n <- 20

test <- numeric(1000)
for (i in 1:1000) {
    x <- runif(n)
    z <- rnorm(n)
    e <- ifelse(runif(n)>p,z, c*z)
    y <- beta*x+e
    fitRQ <- rq(y~x-1)
    w <- x / sum(x^2)
    betatilde <- fitRQ$coef
    s <- sum(abs(fitRQ$res))/n
    test[i] <- betatilde > ( 2 * s * sqrt(sum(w^2)) )
    }

mean(test)

[1] 0.109

p <- 0
c <- 1.5

betavals <- c(0,0.2,0.4,0.6,0.8,1.0)
powersRQ <- numeric(length(betavals))

for (j in (1:length(betavals))){
  beta <- betavals[j]
  for (i in 1:1000) {
    x <- runif(n)
    z <- rnorm(n)
    e <- ifelse(runif(n)>p,z, c*z)
    y <- beta*x+e
    fitRQ <- rq(y~x-1)
    w <- x / sum(x^2)
    betatilde <- fitRQ$coef
    s <- sum(abs(fitRQ$res))/n
    test[i] <- betatilde > 2 * s * sqrt(sum(w^2))
    }
  powersRQ[j] <- mean(test) 
  }

powersRQ

[1] 0.113 0.222 0.343 0.512 0.662 0.783

plot(betavals,powersRQ,ylim=c(0,1))

powersLS <- numeric(length(betavals))

for (j in (1:length(betavals))){
  for (i in 1:1000) {
    beta <- betavals[j]
    x <- runif(n)
    z <- rnorm(n)
    e <- ifelse(runif(n)>p,z, c*z)
    y <- beta*x+e
    fitLS <- lm(y~x-1)
    betahat <- fitLS$coef
    FtestPvalue <- anova(fitLS)$P[1]
    onesidedPvalue <- ifelse(betahat>0, FtestPvalue/2, 1-FtestPvalue/2)
    test[i] <- onesidedPvalue < 0.10
  }
  powersLS[j] <- mean(test) 
}

powersLS

[1] 0.096 0.202 0.389 0.591 0.762 0.875

plot(betavals,powersLS,ylim=c(0,1))

plot(betavals,powersLS,ylim=c(0,1),col="blue")
points(betavals,powersRQ,col="red")