pdf("treasure.pdf") treasure <- read.fwf("C:\\data\\treas1700.dat", width=c(4,11,11,4), col.names=c("year","situados","minted","amerrev")) treasure allmat <- as.matrix(treasure) n <- nrow(allmat) x0 <- matrix(rep(1,n),ncol=1) x <- cbind(x0,allmat[,3]/1000000,allmat[,4]) y <- allmat[,2]/1000000 dfe <- nrow(x)-ncol(x) logy <- log(y) yg <- exp(sum(logy)/n) ipx <- diag(n) - x %*% solve(t(x) %*% x) %*% t(x) lowlam <- -2 hilam <- 2 intrvallam <- 0.1 numlam <- ((hilam-lowlam)/intrvallam)+1 lambda <- matrix(rep(0,numlam),ncol=1) ybc <- matrix(rep(0,n*numlam),ncol=numlam) for (i in 1:numlam) { lambda[i] <- hilam*(i-(numlam+1)/2)/((numlam-1)/2) if (lambda[i]==0) {ybc[,i]<-yg*log(y)} else {ybc[,i]=((y**lambda[i])-1)/(lambda[i]*(yg**(lambda[i]-1)))} } sse <- diag(t(ybc)%*%ipx%*%ybc) lambdasse <- cbind(lambda,sse) for (i in 1:numlam) { if (min(sse)==lambdasse[i,2]) {lam_est <- lambdasse[i,1]} } lam_ci <- min(sse)*(1+(qt(.975,dfe)^2)/dfe) lambdasse print (cbind(lam_est,lam_ci,dfe)) plot(lambda,sse,type='l') abline(h=lam_ci) dev.off()