#****************************************************************************** 
# 1. Spurious Regresson 
# Two random work series regression model
# Simulation based Granger and Newbold (1974)
# "Spurious Regressions in Econometrics 
# 9/26/2017
#------------------------------------------------------------------------------
install.packages("lmtest")
library(lmtest)

# 1) Simple Example 

n = 50 
y = rep(0,n)
x = rep(0,n)
ey = rnorm(n)
ex = rnorm(n)
xrho = 1
yrho = 1
x[1] = 100
y[1] = 100

for (i in 2:n) {
  y[i] = yrho*y[i-1] + ey[i]
  x[i] = xrho*x[i-1] + ex[i]
}
summary(lm( y ~ x))
dwtest(lm( y ~ x))
y <- as.ts(y)
x <- as.ts(x)
ts.plot(y,x,col=c("red", "blue"))



# 2) Simulation Model of Spurious Case

sprious_reg <- function(n,iter) {
  outstat <- matrix(0,iter,6) 
  for (gi in 1:iter){
       y = rep(0,n)
       x = rep(0,n)
      ey = rnorm(n)
      ex = rnorm(n)
      yrho = 1
      xrho = 1
      x[1] = 100
      y[1] = 100
      for (i in 2:n) {
          y[i] = yrho*y[i-1] + ey[i]
          x[i] = xrho*x[i-1] + ex[i]
          }
   out <- summary(lm( y ~ x))
   outstat[gi,1:4] <- out$coefficients[2,]
   outstat[gi,5] <- out$adj.r.squared
   dw <- dwtest(lm( y ~ x))
   outstat[gi,6] <- dw$statistic
   
#   y <- as.ts(y)
#   x <- as.ts(x)
#   ts.plot(y,x,col=c("red", "blue"))
  }
colnames(outstat) <- c("b","se","t","pvalue","adrsq","dwp")
outstat <- as.data.frame(outstat)
return(outstat)  
}

# Customized Tables for interesing statistics 
testresult <- function(dat1){
  table1 <-  summary(dat1$pvalue < .1)
  table2 <-  summary(dat1$pvalue < .05)
  table3 <-  summary(dat1$pvalue < .01)
  table4 <-  summary(dat1$adrsq > .7)
  table5 <-  summary(dat1$dwp < 0.7 )
  table <- rbind(table1,table2,table3,table4,table5)
  rownames(table) <- c("p90","p95","p99","rsq70","dw7")
  return(table)
  }

# Call two function in a line. 
# The first function will generate outstat and the sencond function will create a table
t1 <- Sys.time()
"N = 25 with 100 iteration"
teststat25 <- sprious_reg(25,100000)
testresult(teststat25)
summary(teststat25)

"N = 50 with 100 iteration"
teststat50 <- sprious_reg(50,100000)
testresult(teststat50)
summary(teststat50)

"N = 100 with 100 iteration"
teststat100 <- sprious_reg(100,100000)
testresult(teststat100)
summary(teststat100)

"N = 1000 with 100 iteration"
teststat1k <- sprious_reg(1000,100000)
testresult(teststat1k)
summary(teststat1k)
t2 <- Sys.time()

t2-t1

testresult(teststat25)
testresult(teststat50)
testresult(teststat100)
testresult(teststat1k)