options(object.size=10^100)
##Define FUNCTIONS. We use tukey's tri-weight... no reason... just because
det <- function(M) prod(svd(M, nu = 0, nv = 0)$d)
Window <- function(n, h=1)
{
        aux <- seq(-0.999, 0.999, len = round(n * h))
        Naux1 <- round(n - length(aux))/2
        Naux2 <- n - length(aux) - Naux1
        c(rep(0, Naux1), (1 - abs(aux)^3)^3, rep(0, Naux2))
}
geninv <- function(M){
  aux <- svd(M)
  auxd <- aux$d
  auxd[auxd!=0] <- 1/auxd
  aux$v%*%diag(auxd,nrow=length(auxd))%*%t(aux$u)
}

###READY TO RUN....
sink("tables.tex") ##save ready to latex tables in file
##FIRST CHANGING N and VARIANCE
NN <- 1000 ##Number of replications
Ns <- c(50,100,200,50,100,200)
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
  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 <- 0.5*c(1:K)*log(N) 
  WBICp <- apply(as.array(1:K),1,function(k){
    Y <- as.array(X[,1:k])
    WY <- sweep(Y,1,w,FUN="*")
    0.5*log(det(t(WY)%*%Y)/sigma2) 
  })
  CAICFp <- apply(as.array(1:K),1,function(k){
    Y <- as.array(X[,1:k])
    0.5*(k*log(N) + log(det(t(Y)%*%Y)/N/sigma2) + 2*k )
  })
  CAICFcp <- CAICFp - c(1:K) + (N*(1:K))/(N-(1:K)-2)
  WCAICFp <- WAICp + WBICp 
  res <- array(0,dim=c(NN,8,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)
	wloglik <- 0.5/sigma2*sum(w*residhat^2)
	AIC <- loglik + AICp[k]
	AICc <- loglik + AICcp[k]
	WAIC <- wloglik + WAICp[k]
	BIC <- loglik + BICp[k]
	WBIC <- wloglik + WBICp[k]
	CAICF <- loglik + CAICFp[k]
	CAICFc <- loglik + CAICFcp[k]
	WCAICF <- wloglik + WCAICFp[k]
	c(AIC,AICc,WAIC,BIC,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 <- t(xxx)
  dimnames(xxx) <- list(c("&AIC","&AICc","&WAIC","&BIC","&WBIC","&CAICF","&CAICFc","&WCAICF"),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")}
} 

##NOW KEEP X BETWEEN 0 and 5
BIAS <- 0
Ns <- c(50,100,200,400)
sigma2s <- rep(.25,length(Ns))
for(i in 1:(length(Ns))){
  N <- Ns[i]
##UNFORTUNALTLY I WROTE CODE FOR SIC latter: Here it is:

sink()