ECON546 LAB5 SIMULATION FOR LRT AND WALD TEST

What significance level are we considering in this experiment – 10%, 5% or 1%?

library(mosaic)
#library(forecast)
library(lmtest)
(asymp_crit=qchisq(0.95,1))

[1] 3.841459

#x = seq(1, 5, 0.01 )
library(readr)
xdata <- read_csv("F:/OneDrive/One drive/546-2017/lab/lab5/x.csv")
z = seq(1, 5 , 0.01)
plot( z, dchisq(z, 1),"l")
abline(v = asymp_crit)

eq <- bquote( critical_value == .(asymp_crit))
text(x = asymp_crit+0.6, y = 0.05, eq)

In which line of the code do we assign the null hypothesis to be true?

rho=0.0

What sample size is being considered?

n=10

Which line(s) would we alter if we wanted to consider a sample size of (say) n = 200?

n = 200

the percentage rejection rates for the case where THE NULL HYPOTHESIS IS TRUE - so these are the true sizes (significance levels) of the tests,

set.seed(654321)
nrep=2000
n=10
c1=1 
c2=4
rho=0.0
#x=rnorm(n)+(1:n)/n
x = xdata$x[1:n]
#str(xdata$x)
#y = numeric(20000) 
 
lrtv = numeric(nrep)
waldv = numeric(nrep)
power_lrt=0.0 
power_wald=0.0 
asymp_crit=qchisq(0.95,1)
#'============================='
for (i in 1: nrep ){
    
  e = numeric(n)    
  #genr e=rho*e(-1)+@rnorm 
  #e = arima.sim(model=list(ar=rho),n=n) 
  
  # or
  
  for(j in 2:n){
    e[j] = rho*e[j-1] + rnorm(1) 
  }
  #arima(e, order = c(1, 0, 0)) 
  #acf(e, 10)
  #pacf(e, 10)
  y=c1+c2*x+e 
  
  
  #equation eq1.ls y c x 
  fit1 <- lm(y~ x)
  #summary(fit1)  
  
  #loglr=@logl 
  
  loglr=logLik(fit1)  
  
  #http://stackoverflow.com/questions/7233288/non-stationary-seasonal-ar-part-from-css-error-in-r
  #equation eq2.ls(arma=ml) y c x ar(1) 
  fit2 <- arima(y, xreg=x, order=c(1,0,0),method="ML")
  #summary(fit2)
  #str(coeftest(fit2) )
  
  #scalar loglu=@logl 
  
  loglu=logLik(fit2)
  
  lrtv[i]=-2*(loglr-loglu)
  
  
  #'if lrtv(!i)>asymp_crit then
  if (lrtv[i]>asymp_crit){
  
  
  #if (lrtv[i]>crit_lrt){
    power_lrt=power_lrt+1 
    
  } 
   
  
  
  waldv[i]=coeftest(fit2)[1,3]^2 
  
  
  #'if waldv(!i)>asymp_crit
  
  if (waldv[i]>asymp_crit){
 
  
  #if (waldv[i]>crit_lrt){
    
    power_wald=power_wald+1
    
    } 
}
#'============================='
(power_lrt=power_lrt/nrep*100)

[1] 7.6

(power_wald=power_wald/nrep*100)

[1] 13.3

(crit_lrt=quantile(lrtv,0.95) )

     95% 
4.815507 

hist(lrtv)

(crit_wald=quantile(waldv,0.95))

     95% 
7.825018 

hist(waldv)

Make a function

sizeAdjust = function(  nseed=654321, rho=0.0, nrep = 2000, n = 10, alpha = 0.05){
  
    set.seed(nseed)
    
    #nrep=nrep
    #n=10
    
    c1=1 
    c2=4
    #rho=0.0
    #x=rnorm(n)+(1:n)/n
    x = xdata$x[1:n]
    
    #y = numeric(20000) 
     
    lrtv = numeric(nrep)
    waldv = numeric(nrep)
    
    power_lrt=0.0 
    
    power_wald=0.0 
    asymp_crit=qchisq(1-alpha,1)
    
    #'============================='
    
    
    
    for (i in 1: nrep ){
    
        
      e = numeric(n)    
      #genr e=rho*e(-1)+@rnorm 
      #e = arima.sim(model=list(ar=rho),n=n) 
      
      # or
      
      for(j in 2:n){
        e[j] = rho*e[j-1] + rnorm(1) 
      }
      #arima(e, order = c(1, 0, 0)) 
      #acf(e, 10)
      #pacf(e, 10)
      y=c1+c2*x+e 
      
      
      #equation eq1.ls y c x 
      fit1 <- lm(y~ x)
      #summary(fit1)  
      
      #loglr=@logl 
      
      loglr=logLik(fit1)  
      
      #http://stackoverflow.com/questions/7233288/non-stationary-seasonal-ar-part-from-css-error-in-r
      #equation eq2.ls(arma=ml) y c x ar(1) 
      fit2 <- arima(y, xreg=x, order=c(1,0,0),method="ML")
      #summary(fit2)
      #str(coeftest(fit2) )
      
      #scalar loglu=@logl 
      
      loglu=logLik(fit2)
      
      lrtv[i]=-2*(loglr-loglu)
      
      
      #'if lrtv(!i)>asymp_crit then
      if (lrtv[i]>asymp_crit){
      
      
      #if (lrtv[i]>crit_lrt){
        power_lrt=power_lrt+1 
        
      } 
       
      
      
      waldv[i]=coeftest(fit2)[1,3]^2 
      
      
      #'if waldv(!i)>asymp_crit
      
      if (waldv[i]>asymp_crit){
     
      
      #if (waldv[i]>crit_lrt){
        
        power_wald=power_wald+1
        
        } 
    
    }
       #'============================='
    
    power_lrt=power_lrt/nrep*100
    
    
    power_wald=power_wald/nrep*100
    
    crit_lrt=quantile(lrtv,0.95) 
    crit_wald=quantile(waldv,0.95)
    
    
    return(c(power_lrt, power_wald,crit_lrt,crit_wald ))
}

Computer power

computePower =function( nseed=654321, rho=0.0, nrep = 2000, n = 10, alpha = 0.05， crit_lrt = 4.815507， crit_wald = 7.825018){
  
    set.seed(nseed)
    
    #nrep=nrep
    #n=10
    
    c1=1 
    c2=4
    #rho=0.0
    #x=rnorm(n)+(1:n)/n
    x = xdata$x[1:n]
    
    #y = numeric(20000) 
     
    lrtv = numeric(nrep)
    waldv = numeric(nrep)
    
    power_lrt=0.0 
    
    power_wald=0.0 
    asymp_crit=qchisq(1-alpha,1)
    
    #'============================='
    
    
    
    for (i in 1: nrep ){
    
        
      e = numeric(n)    
      #genr e=rho*e(-1)+@rnorm 
      #e = arima.sim(model=list(ar=rho),n=n) 
      
      # or
      
      for(j in 2:n){
        e[j] = rho*e[j-1] + rnorm(1) 
      }
      #arima(e, order = c(1, 0, 0)) 
      #acf(e, 10)
      #pacf(e, 10)
      y=c1+c2*x+e 
      
      
      #equation eq1.ls y c x 
      fit1 <- lm(y~ x)
      #summary(fit1)  
      
      #loglr=@logl 
      
      loglr=logLik(fit1)  
      
      #http://stackoverflow.com/questions/7233288/non-stationary-seasonal-ar-part-from-css-error-in-r
      #equation eq2.ls(arma=ml) y c x ar(1) 
      fit2 <- arima(y, xreg=x, order=c(1,0,0),method="ML")
      #summary(fit2)
      #str(coeftest(fit2) )
      
      #scalar loglu=@logl 
      
      loglu=logLik(fit2)
      
      lrtv[i]=-2*(loglr-loglu)
      
      
      #'if lrtv(!i)>asymp_crit then
      if (lrtv[i]>crit_lrt){
      
      
      #if (lrtv[i]>crit_lrt){
        power_lrt=power_lrt+1 
        
      } 
       
      
      
      waldv[i]=coeftest(fit2)[1,3]^2 
      
      
      #'if waldv(!i)>asymp_crit
      
      if (waldv[i]>asymp_crit){
     
      
      #if (waldv[i]>crit_lrt){
        
        power_wald=power_wald+1
        
        } 
    
    }
       #'============================='
    
    power_lrt=power_lrt/nrep*100
    
    
    power_wald=power_wald/nrep*100
    
    #crit_lrt=quantile(lrtv,0.95) 
    #crit_wald=quantile(waldv,0.95)
    
    
    return(c(power_lrt, power_wald))
}

Plot power curve

sizeAdjust(nseed=654321, rho=0.0, nrep = 2000, n = 10, alpha = 0.05)

                          95%       95% 
 7.600000 13.300000  4.815507  7.825018 

Change parameters

sizeAdjust(nseed=654321, rho=0.0, nrep = 2000, n = 200, alpha = 0.05)

                       95%      95% 
5.650000 5.700000 4.034827 4.077705 

Change parameters

sizeAdjust(nseed=654321, rho=0.0, nrep = 2000, n = 500, alpha = 0.05)

                       95%      95% 
5.350000 5.400000 3.913753 3.930252 

Change parameters

sizeAdjust(nseed=654321, rho=0.0, nrep = 2000, n = 1000, alpha = 0.05)

                       95%      95% 
4.550000 4.550000 3.635447 3.641662 

Change parameters

computePower(nseed=654321, rho=-0.1,  nrep = 2000, n = 10, alpha = 0.05)

[1]  7.05 16.40

Change parameters

computePower(nseed=654321, rho=-0.3,  nrep = 2000, n = 10, alpha = 0.05)

[1] 15.05 28.75

Change parameters

computePower(nseed=654321, rho=-0.5,  nrep = 2000, n = 10, alpha = 0.05)

[1] 31.75 49.20

Change parameters

computePower(nseed=654321, rho=-0.7,  nrep = 2000, n = 10, alpha = 0.05)

[1] 53.60 69.55

Change parameters

computePower(nseed=654321, rho=0.1,  nrep = 2000, n = 10, alpha = 0.05)

[1]  3.7 11.4

Change parameters

computePower(nseed=654321, rho=0.3,  nrep = 2000, n = 10, alpha = 0.05)

[1]  3.85 12.35

Change parameters

computePower(nseed=654321, rho=0.5,  nrep = 2000, n = 10, alpha = 0.05)

[1]  6.4 18.8

Change parameters

computePower(nseed=654321, rho=0.7,  nrep = 2000, n = 10, alpha = 0.05)

[1] 11.75 29.55