######## THIS IS THE SIMULATION PRESENTED IN SECTION 5.3 ####### IT is the LOCAL REGRESSION SIMULATION ## THE SICf is at the end. options(object.size=10^100) ##FUNCTIONS: det returns the determinant of matrix, geninv is generlized ##and Window defines the window function... we use gaussian here for ease ##of computation det <- function(M) prod(svd(M, nu = 0, nv = 0)$d) geninv <- function(M){ aux <- svd(M) auxd <- aux$d auxd[auxd!=0] <- 1/auxd aux$v%*%diag(auxd,nrow=length(auxd))%*%t(aux$u) } Window <- function(o,middle=pi/2,sigma=5){ dnorm((o-middle)/max(abs(o-middle))*sigma) } ###READY TO RUN.... Save ready for latex table to table.tex sink("tables.tex") ##FIRST CHANGING N and VARIANCE NN <- 1000##Number of replications cat(rep("#",20),"\n","SIM 2\n",rep("#",20),"\n") Ns <- rep(c(50,100,200)+1,2) #different values of n sigma2s <- c(rep(1,length(Ns)/2),rep(5,length(Ns)/2))#different values of sigma for(i in 1:(length(Ns))){ N <- Ns[i] sigma2 <- sigma2s[i] cat(N,",",sigma2,"\n") ##we leave a hole out so that constant and line dont give the same estimate HALF <- (N-1)/2; HOLE <- round(.06*HALF) o1 <- c(runif(HALF-HOLE,0,pi/2-.2),runif(HOLE,pi/2-.1,pi/2)) o <- sort(c(o1,pi/2,runif(HALF,pi/2,3/2*pi))) w <- Window(o) w <- w/sum(w)*N W0 <- sum(w) truemean <- 2*cos(o) ##THIS IS THE TRUE MEAN s(x) o <- o*10 #make o bigger to avoid problems with geninv K <- 4 X <- apply(as.array(0:(K-1)),1,function(k) (o^k)) AICp <- c(1:K) AICcp <- AICp + (c(1:K) + 1)*(c(1:K) + 2)/(N - c(1:K) - 2) WAICp <- apply(as.array(1:K),1,function(k){ Y <- as.array(X[,1:k]) WY <- sweep(Y,1,w,FUN="*") aux <- geninv(t(Y)%*%WY) W2Y <- sweep(Y,1,w^2,FUN="*") sum(diag(aux%*%(t(W2Y)%*%Y))) }) BICp <- c(1:K)*log(N) SICp <- apply(as.array(1:K),1,function(k){ Y <- as.array(X[,1:k]) log(det(t(Y)%*%Y)) }) WBICp <- apply(as.array(1:K),1,function(k){ Y <- as.array(X[,1:k]) WY <- sweep(Y,1,w,FUN="*") log(det(t(WY)%*%Y)) }) CAICFp <- apply(as.array(1:K),1,function(k){ Y <- as.array(X[,1:k]) (log(det(t(Y)%*%Y)) + 2*k) }) CAICFcp <- CAICFp - c(1:K) + (N*(1:K))/(N-(1:K)-2) WCAICFp <- WAICp + WBICp res <- array(0,dim=c(NN,10,K)) for(j in 1:NN) { y <- truemean + rnorm(N,0,sqrt(sigma2)) res[j,,] <- apply(as.array(1:K),1,function(k){ auxfit <- lm(y~X[,1:k]-1,weight=w) residhat <- auxfit$resid mse <- (auxfit$fitted[(N+1)/2] - truemean[(N+1)/2])^2 sigma2hat <- sum(residhat^2)/(N-AICp[k]) wsigma2hat <- sum(w*residhat^2)/(W0-WAICp[k]) AIC <- N*log(sigma2hat) + AICp[k] AICc <- N*log(sigma2hat) + AICcp[k] WAIC <- W0*log(wsigma2hat) + WAICp[k] BIC <- (N-1)*log(sigma2hat) + BICp[k] SIC <- (N-1)*log(sigma2hat) + SICp[k] WBIC <- (W0-1)*log(wsigma2hat) + WBICp[k] CAICF <- (N-1)*log(sigma2hat) + CAICFp[k] CAICFc <- (N-1)*log(sigma2hat) + CAICFcp[k] WCAICF <- (W0-1)*log(wsigma2hat) + WCAICFp[k] c(mse,AIC,AICc,WAIC,BIC,SIC,WBIC,CAICF,CAICFc,WCAICF) }) } xx <- apply(res,c(1,2),function(x) order(x)[1]) xxx <- apply(xx,2,function(x) table(factor(x,levels=c(1:K)))/NN*100) xxx[,1] <- signif(apply(res[,1,],2,mean),8) xxx <- t(xxx) dimnames(xxx) <- list(c("&MSE","&AIC","&AICc","&WAIC","&BIC","&SICf","&WBIC","&CAICF","&CAICFc","&WCAICF"),as.character(c(1:K))) cat("\n") for(h in 1:dim(xxx)[1]){ ##HERE WE PRINT OUT THE RESULT if(h==5) cat("$n=",N,"$") if(h==6) cat("$\sigma^2=",sigma2,"$") cat(dimnames(xxx)[[1]][h]) for(g in 1:dim(xxx)[2]) cat("&",xxx[h,g]) cat("&\n")} } sink() sink("SICf.tex") ##FIRST CHANGING N and VARIANCE #BIAS <- 6/24 #BIAS <- 25 NN <- 1000 ##Number of replications Ns <- c(50,100,200,50,100,200) #sigmas2s <- rep(.25,3) sigma2s <- c(.25,.5,1,5,5,5) for(i in 1:(length(Ns))){ N <- Ns[i] cat(N,"\n") K <- 7 trueK <- 4 o <- c(0:(N-1)) N <- length(o) beta <- matrix(c(1,5,-1.25,0.15,0,0,0),ncol=1) X <- apply(as.array(0:(K-1)),1,function(k) (o^k)) sigma2 <- sigma2s[i] truemean <- (X%*%beta) w <- Window(N) w <- w/sum(w)*N BICp <- apply(as.array(1:K),1,function(k){ Y <- as.array(X[,1:k]) 0.5*log(det(t(Y)%*%Y)/sigma2) }) res <- array(0,dim=c(NN,1,K)) for(j in 1:NN) { y <- truemean + rnorm(N,0,sqrt(sigma2)) res[j,,] <- apply(as.array(1:K),1,function(k){ residhat <- lm(y~X[,1:k]-1,weight=w)$resid loglik <- 0.5/sigma2*sum(residhat^2) BIC <- loglik + BICp[k] BIC }) } xx <- apply(res,c(1,2),function(x) order(x)[1]) xxx <- apply(xx,2,function(x) table(factor(x,levels=c(1:K)))/NN*100) xxx <- t(xxx) dimnames(xxx) <- list(c("&SICf"),as.character(c(1:K))) for(h in 1:dim(xxx)[1]){ if(h==4) cat("$n=",N,"$") if(h==5) cat("$\sigma^2=",sigma2,"$") cat(dimnames(xxx)[[1]][h]) for(g in 1:dim(xxx)[2]) cat("&",xxx[h,g]) cat("&\n")} } sink()