# R 3.0.1
#  Ex4_SpurReg_R.txt
## Monte Carlo simulation for the 
## spurious regression problem

rm(list=ls())
library("tseries")  
library("zoo")
library("dynlm")

nsim <- 100                                # no of replications/simulations
# nsim <- 1
n <- 100                                   # choose sample size

b1 <- r21 <- tstat1 <- c(1:nsim)*0;    # initialization
b0 <- r20 <- tstat0 <- c(1:nsim)*0;    # initialization

# simulation loop
for (isim in 1:nsim) 
  { 
 dx <- rnorm(n); dy <- rnorm(n);
 x  <- dx;           y <- dy;
 for (t in 2:n) { 
    x[t]  <- x[t-1] + dx[t]; y[t]  <- y[t-1] + dy[t];
                   }
  xy   <- cbind(x,y);              # x,   y two independent random walks
  dxy <- cbind(dx,dy);             # dx, dy two independent white noises
# plot only for nsim = 1
 if ( nsim == 1) {
    plot(dx,dy) ; plot(dxy);
    plot(xy);       plot(x,y)
                      }
#

# estimate LS
 mod0 <- lm(dy ~ dx + 1)
 # summary(mod0)                      # see       help(summary.lm)
 b0[isim]        <- summary(mod0)$coef[2,1] # slope
 tstat0[isim]   <- summary(mod0)$coef[2,3] # t-statistic of slope
 r20[isim]       <- summary(mod0)$r.squared 
#
 mod1 <- lm(y ~ x + 1)
 # summary(mod1)                      # see       help(summary.lm)
 b1[isim]        <- summary(mod1)$coef[2,1] # slope
 tstat1[isim]   <- summary(mod1)$coef[2,3] # t-statistic of slope
 r21[isim]       <- summary(mod1)$r.squared 
  }
# end simulation loop

plot(density(b0),      type="h")
plot(density(tstat0), type="h")
plot(density(r20),     type="h")

plot(ecdf( abs(tstat0) ))       # empirical distribution function
                                             # check proportion > 1.96 
#
plot(density(b1),      type="h")
plot(density(tstat1), type="h")
plot(density(r21),     type="h")

plot(ecdf( abs(tstat1) ))       # empirical distribution function
                                             # check proportion > 1.96 


## EXERCISE: run this script 
# (1) as it is
# (2) repeatedly with nsim <- 1 in order to inspect 
#      the plots of x and y, dx, dy
#      the scatterplot (x,y), scatterplot(dx,dy)
#      single estimation results (R2, t-value, DW)
# (3) repeat with n <- 400, 1600, ... 
# 
# Compare the densities of b, tstat, r2 and abs(tstat) with their 
# theoretical properties.