####################################################################

##Execute following on dos command prompt (omit the "##"):
##"S:\\R\\rw2001\\bin\\Rterm.exe" -q --no-restore --no-save < coint_test.r > R.out

####################################################################

options(echo = FALSE)
library(fSeries);
library(MASS);
suppressWarnings(library(urca));

png(file="plot%03d.png", width = 680, height = 450, pointsize = 12, bg="transparent", res = 200)

#--------------------------------------------------------------------
# Create 2 series: "x", "y" starting Jan 1960 with random values assigned
# for 100 years from start for each. 12*100 = 1200 values. Store in time
# series object "z".

numyears <- 100;
startyr <- 1960;
startper <- 1;
freq <- 12;

# Simulated cointegrated series pair:
# random walk = I(1) process:
xc <- cumsum(rnorm(numyears*freq));
# linear combination of the above, also = I(1):
yc <- 20.0 + 10.0 * xc + rnorm(length(xc), mean=0, sd=5);
zc <- ts( cbind(xc,yc), start=c(startyr, startper), frequency=freq );

# Simulated non-cointegrated series pair:
# stationary pure (pseudo) white noise = I(0):
xnc <- rnorm(numyears*freq);
# random walk = nonstationary I(1) process:
ync <- cumsum(rnorm(numyears*freq));
znc <- ts( cbind(xnc, ync), start=c(startyr, startper), frequency=freq );

All_Series <- ts( cbind(xc,yc,xnc,ync), start=c(startyr, startper), frequency=freq);
plot(All_Series);

# Compute first differences of above series:
xc_diff <- diff(xc);
yc_diff <- diff(yc);
xnc_diff <- diff(xnc);
ync_diff <- diff(ync);

#--------------------------------------------------------------------
# Unit root test 1: Augmented Dickey-Fuller test on all simulated series.
# Ho: non-stationary. Perform on input series and on first diff. series.
# Fix maximum lag in regression to 2 (may not be optimal!)

maxlag <- 2;

adf_xc <- suppressWarnings(tsadfTest(xc, alternative = "stationary", k=maxlag));
adf_yc <- suppressWarnings(tsadfTest(yc, alternative = "stationary", k=maxlag));
adf_xnc <- suppressWarnings(tsadfTest(xnc, alternative = "stationary", k=maxlag));
adf_ync <- suppressWarnings(tsadfTest(ync, alternative = "stationary", k=maxlag));

adf_stat_xc <- adf_xc@test$statistic;
adf_stat_yc <- adf_yc@test$statistic;
adf_stat_xnc <- adf_xnc@test$statistic;
adf_stat_ync <- adf_ync@test$statistic;
adf_pvalue_xc <- adf_xc@test$p.value;
adf_pvalue_yc <- adf_yc@test$p.value;
adf_pvalue_xnc <- adf_xnc@test$p.value;
adf_pvalue_ync <- adf_ync@test$p.value;

title <- "Augmented Dickey-Fuller Unit Root Test on input series:";
write(paste("\n*",title), file="");
write(cat("* ",(rep("-",trunc(0.5*(nchar(title) + 2)))), sep="-"), file="");
write(paste("* Series xc: ADF_STATISTIC =",adf_stat_xc), file="");
write(paste("* Series xc: P-VALUE (H0: non-stationarity) =",adf_pvalue_xc), file="");
write(paste("* Series yc: ADF_STATISTIC =",adf_stat_yc), file="");
write(paste("* Series yc: P-VALUE (H0: non-stationarity) =",adf_pvalue_yc), file="");
write(paste("* Series xnc: ADF_STATISTIC =",adf_stat_xnc), file="");
write(paste("* Series xnc: P-VALUE (H0: non-stationarity) =",adf_pvalue_xnc), file="");
write(paste("* Series ync: ADF_STATISTIC =",adf_stat_ync), file="");
write(paste("* Series ync: P-VALUE (H0: non-stationarity) =",adf_pvalue_ync), file="");

# Run on first differences:

adf_xc_diff <- suppressWarnings(tsadfTest(xc_diff, alternative = "stationary", k=maxlag));
adf_yc_diff <- suppressWarnings(tsadfTest(yc_diff, alternative = "stationary", k=maxlag));
adf_xnc_diff <- suppressWarnings(tsadfTest(xnc_diff, alternative = "stationary", k=maxlag));
adf_ync_diff <- suppressWarnings(tsadfTest(ync_diff, alternative = "stationary", k=maxlag));

adf_stat_xc_diff <- adf_xc_diff@test$statistic;
adf_stat_yc_diff <- adf_yc_diff@test$statistic;
adf_stat_xnc_diff <- adf_xnc_diff@test$statistic;
adf_stat_ync_diff <- adf_ync_diff@test$statistic;
adf_pvalue_xc_diff <- adf_xc_diff@test$p.value;
adf_pvalue_yc_diff <- adf_yc_diff@test$p.value;
adf_pvalue_xnc_diff <- adf_xnc_diff@test$p.value;
adf_pvalue_ync_diff <- adf_ync_diff@test$p.value;

title <- "Augmented Dickey-Fuller Unit Root Test on first-differenced series:";
write(paste("\n*",title), file="");
write(cat("* ",(rep("-",trunc(0.5*(nchar(title) + 2)))), sep="-"), file="");
write(paste("* Series xc_diff: ADF_STATISTIC =",adf_stat_xc_diff), file="");
write(paste("* Series xc_diff: P-VALUE (H0: non-stationarity) =",adf_pvalue_xc_diff), file="");
write(paste("* Series yc_diff: ADF_STATISTIC =",adf_stat_yc_diff), file="");
write(paste("* Series yc_diff: P-VALUE (H0: non-stationarity) =",adf_pvalue_yc_diff), file="");
write(paste("* Series xnc_diff: ADF_STATISTIC =",adf_stat_xnc_diff), file="");
write(paste("* Series xnc_diff: P-VALUE (H0: non-stationarity) =",adf_pvalue_xnc_diff), file="");
write(paste("* Series ync_diff: ADF_STATISTIC =",adf_stat_ync_diff), file="");
write(paste("* Series ync_diff: P-VALUE (H0: non-stationarity) =",adf_pvalue_ync_diff), file="");

#--------------------------------------------------------------------
# Unit root test 2: Phillips-Perron test on all simulated series.
# Ho: non-stationary. Perform on input series and on first diff. series.

pp_xc <- suppressWarnings(tsppTest(xc, alternative="stationary", type="Z(alpha)", lshort=TRUE));
pp_yc <- suppressWarnings(tsppTest(yc, alternative="stationary", type="Z(alpha)", lshort=TRUE));
pp_xnc <- suppressWarnings(tsppTest(xnc, alternative="stationary", type="Z(alpha)", lshort=TRUE));
pp_ync <- suppressWarnings(tsppTest(ync, alternative="stationary", type="Z(alpha)", lshort=TRUE));

pp_stat_xc <- pp_xc@test$statistic;
pp_stat_yc <- pp_yc@test$statistic;
pp_stat_xnc <- pp_xnc@test$statistic;
pp_stat_ync <- pp_ync@test$statistic;
pp_pvalue_xc <- pp_xc@test$p.value;
pp_pvalue_yc <- pp_yc@test$p.value;
pp_pvalue_xnc <- pp_xnc@test$p.value;
pp_pvalue_ync <- pp_ync@test$p.value;

title <- "Phillips-Perron Unit Root Test on input series:";
write(paste("\n*",title), file="");
write(cat("* ",(rep("-",trunc(0.5*(nchar(title) + 2)))), sep="-"), file="");
write(paste("* Series xc: Z_STATISTIC =",pp_stat_xc), file="");
write(paste("* Series xc: P-VALUE (H0: non-stationarity) =",pp_pvalue_xc), file="");
write(paste("* Series yc: Z_STATISTIC =",pp_stat_yc), file="");
write(paste("* Series yc: P-VALUE (H0: non-stationarity) =",pp_pvalue_yc), file="");
write(paste("* Series xnc: Z_STATISTIC =",pp_stat_xnc), file="");
write(paste("* Series xnc: P-VALUE (H0: non-stationarity) =",pp_pvalue_xnc), file="");
write(paste("* Series ync: Z_STATISTIC =",pp_stat_ync), file="");
write(paste("* Series ync: P-VALUE (H0: non-stationarity) =",pp_pvalue_ync), file="");

# Run on first differences:

pp_xc_diff <- suppressWarnings(tsppTest(xc_diff, alternative="stationary", type="Z(alpha)", lshort=TRUE));
pp_yc_diff <- suppressWarnings(tsppTest(yc_diff, alternative="stationary", type="Z(alpha)", lshort=TRUE));
pp_xnc_diff <- suppressWarnings(tsppTest(xnc_diff, alternative="stationary", type="Z(alpha)", lshort=TRUE));
pp_ync_diff <- suppressWarnings(tsppTest(ync_diff, alternative="stationary", type="Z(alpha)", lshort=TRUE));

pp_stat_xc_diff <- pp_xc_diff@test$statistic;
pp_stat_yc_diff <- pp_yc_diff@test$statistic;
pp_stat_xnc_diff <- pp_xnc_diff@test$statistic;
pp_stat_ync_diff <- pp_ync_diff@test$statistic;
pp_pvalue_xc_diff <- pp_xc_diff@test$p.value;
pp_pvalue_yc_diff <- pp_yc_diff@test$p.value;
pp_pvalue_xnc_diff <- pp_xnc_diff@test$p.value;
pp_pvalue_ync_diff <- pp_ync_diff@test$p.value;

title <- "Phillips-Perron Unit Root Test on first-differenced series:";
write(paste("\n*",title), file="");
write(cat("* ",(rep("-",trunc(0.5*(nchar(title) + 2)))), sep="-"), file="");
write(paste("* Series xc_diff: Z_STATISTIC =",pp_stat_xc_diff), file="");
write(paste("* Series xc_diff: P-VALUE (H0: non-stationarity) =",pp_pvalue_xc_diff), file="");
write(paste("* Series yc_diff: Z_STATISTIC =",pp_stat_yc_diff), file="");
write(paste("* Series yc_diff: P-VALUE (H0: non-stationarity) =",pp_pvalue_yc_diff), file="");
write(paste("* Series xnc_diff: Z_STATISTIC =",pp_stat_xnc_diff), file="");
write(paste("* Series xnc_diff: P-VALUE (H0: non-stationarity) =",pp_pvalue_xnc_diff), file="");
write(paste("* Series ync_diff: Z_STATISTIC =",pp_stat_ync_diff), file="");
write(paste("* Series ync_diff: P-VALUE (H0: non-stationarity) =",pp_pvalue_ync_diff), file="");

rm(xc_diff, yc_diff, xnc_diff, ync_diff);

#--------------------------------------------------------------------
# Unit roots conclusions:

integ_flag_xc <- 0;
integ_flag_yc <- 0;
integ_flag_xnc <- 0;
integ_flag_ync <- 0;

# series xc:
if( (adf_pvalue_xc > 0.05) && (pp_pvalue_xc > 0.05) &&
    (adf_pvalue_xc_diff < 0.05) && (pp_pvalue_xc_diff < 0.05) ) {
  write(paste("\n*** Series xc = integrated I(1) process [at >5% level; see pvalues above]"),file="");
  integ_flag_xc <- 1;  
} else {
  write(paste("\n*** Series xc does not appear to be an integrated I(1) process [at <5% level; see pvalues above]"),file="");
}

# series yc:
if( (adf_pvalue_yc > 0.05) && (pp_pvalue_yc > 0.05) &&
    (adf_pvalue_yc_diff < 0.05) && (pp_pvalue_yc_diff < 0.05) ) {
  write(paste("*** Series yc = integrated I(1) process [at >5% level; see pvalues above]"),file="");
  integ_flag_yc <- 1;  
} else {
  write(paste("*** Series yc does not appear to be integrated I(1) process [at <5% level; see pvalues above]"),file="");
}

# series xnc:
if( (adf_pvalue_xnc > 0.05) && (pp_pvalue_xnc > 0.05) &&
    (adf_pvalue_xnc_diff < 0.05) && (pp_pvalue_xnc_diff < 0.05) ) {
  write(paste("\n*** Series xnc = integrated I(1) process [at >5% level; see pvalues above]"),file="");
  integ_flag_xnc <- 1;  
} else {
  write(paste("*** Series xnc does not appear to be an integrated I(1) process [at <5% level; see pvalues above]"),file="");
}

# series ync:
if( (adf_pvalue_ync > 0.05) && (pp_pvalue_ync > 0.05) &&
    (adf_pvalue_ync_diff < 0.05) && (pp_pvalue_ync_diff < 0.05) ) {
  write(paste("*** Series ync = integrated I(1) process [at >5% level; see pvalues above]"),file="");
  integ_flag_ync <- 1;  
} else {
  write(paste("*** Series ync does not appear to be integrated I(1) process [at <5% level; see pvalues above]"),file="");
}

#--------------------------------------------------------------------
# Regress series Yt on series Xt if both are I(1): Yt = B.Xt + A.
# Use "robust" linear regression using "rlm" function (outlier resistant).
# Try "MM" method first, if fails (from size of returned object
# containing just error message string), re-execute with "M" method.
# Otherwise let fail.

# Ensure first that series have same length!

if( integ_flag_xc && integ_flag_yc ) {
  lin_co <- try(rlm(yc ~ xc + 1, method="MM"), silent=TRUE);
  if( length(lin_co) == 1 ) {
    lin_co <- rlm(yc ~ xc + 1, method="M");
  }
  resid_co <- residuals(lin_co);
}

if( integ_flag_xnc && integ_flag_ync ) {
  lin_noco <- try(rlm(ync ~ xnc + 1, method="MM"), silent=TRUE);
  if( length(lin_noco) == 1 ) {
    lin_noco <- rlm(ync ~ xnc + 1, method="M");
  }
  resid_noco <- residuals(lin_noco);
}

#--------------------------------------------------------------------
# Cointegration test 1: Augmented Engle-Granger (AEG) test on residuals
# of cointegrating regression: Ut = Yt - B.Xt - A. Use critical tau statistics
# for cointegration tests from MacKinnon, J.G. 1996.

N <- numyears * freq;

critval_0.01 <- -4.3266 + (-15.531/N) + (-34.03*(1/N)**2);
critval_0.05 <- -3.7809 + (-9.421/N) + (-15.06*(1/N)**2);
critval_0.10 <- -3.4959 + (-7.203/N) + (-4.01*(1/N)**2);

title <- "Cointegration Test 1: Augmented Engle-Granger (AEG) test on residuals of regression: Y ~ X:";
write(paste("\n*",title), file="");
write(cat("* ",(rep("-",trunc(0.5*(nchar(title) + 2)))), sep="-"), file="");

if( integ_flag_xc && integ_flag_yc ) {
  adf_resid_coint <- suppressWarnings(tsadfTest(resid_co, alternative="stationary", k=maxlag));
  adf_stat <- adf_resid_coint@test$statistic;

  write(paste("* Residuals series [u = yc - B.xc - A]  ADF_STATISTIC =",adf_stat), file="");
  write(paste("* Critical values [p%] for AEG test =",round(critval_0.01,3),"[1%]",round(critval_0.05,3),"[5%]",round(critval_0.10,3),"[10%]"), file="");

  if( abs(adf_stat) > abs(critval_0.01) ) {
    write(paste("\n*** Series pair (xc,yc) is cointegrating [H0: no-cointegration rejected at <1% level]"),file="");
  } else if( abs(adf_stat) > abs(critval_0.05) ) {
    write(paste("\n*** Series pair (xc,yc) is cointegrating [H0: no-cointegration rejected at <5% level]"),file="");
  } else {
    write(paste("\n*** Series pair (xc,yc) is not cointegrating [H0: no-cointegration accepted at >5% level]"),file="");
  }

} else {
  write(paste("\n*** Series pair (xc,yc) is not cointegrating since series are not integrated with same order."),file="");
}

if( integ_flag_xnc && integ_flag_ync ) {
  adf_resid_nocoint <- suppressWarnings(tsadfTest(resid_noco, alternative="stationary", k=maxlag));
  adf_stat <- adf_resid_nocoint@test$statistic;

  write(paste("* Residuals series [u = ync - B.xnc - A]  ADF_STATISTIC =",adf_stat), file="");
  write(paste("* Critical values [p%] for AEG test =",round(critval_0.01,3),"[1%]",round(critval_0.05,3),"[5%]",round(critval_0.10,3),"[10%]"), file="");

  if( abs(adf_stat) > abs(critval_0.01) ) {
    write(paste("\n*** Series pair (xnc,ync) is cointegrating [H0: no-cointegration rejected at <1% level]"),file="");
  } else if( abs(adf_stat) > abs(critval_0.05) ) {
    write(paste("\n*** Series pair (xnc,ync) is cointegrating [H0: no-cointegration rejected at <5% level]"),file="");
  } else {
    write(paste("\n*** Series pair (xnc,ync) is not cointegrating [H0: no-cointegration accepted at >5% level]"),file="");
  }

} else {
  write(paste("*** Series pair (xnc,ync) is not cointegrating since series are not integrated with same order."),file="");
}

#--------------------------------------------------------------------
# Cointegration test 2: Johansen's VAR-ECM procedure for
# significantly non-zero eigenvalues in long run matrix.

title <- "Cointegration Test 2: Johansen's VAR procedure for significantly non-zero eigenvalues:";
write(paste("\n*",title), file="");
write(cat("* ",(rep("-",trunc(0.5*(nchar(title) + 2)))), sep="-"), file="");

#---- xc, yc (cointegrated) series pair: ----

if( integ_flag_xc && integ_flag_yc ) {
  # Trace test:
  vecm_trace_coint <- ca.jo(zc, type="trace", constant=FALSE, K=2, spec="longrun", season=NULL, dumvar=NULL, ctable="A2");
  eigenvalues <- vecm_trace_coint@lambda;
  tr_statistics <- vecm_trace_coint@teststat;
  tr_cvals <- vecm_trace_coint@cval;

  write(paste("* With max. lag = 2 and linear trend included."), file="");
  write(paste("* Eigenvalue (e1) =", round(eigenvalues[[2]],8),";","Eigenvalue (e2) =",round(eigenvalues[[1]],8)), file="");
  write(paste("\n* Trace statistic[e1 + e2] (H0: r=0 cointegrating vectors vs. H1: r=2) =", round(tr_statistics[[2]],4)), file="");
  write(paste("* Trace statistic[e1] (H0: r<=1 cointegrating vector vs. H1: r=2) =", round(tr_statistics[[1]],4)), file="");
  write(paste("* Critical trace test values [r=0]: ", tr_cvals[[6]],"[1%]",tr_cvals[[4]],"[5%]",tr_cvals[[2]],"[10%]"), file="");
  write(paste("* Critical trace test values [r<=1]: ", tr_cvals[[5]],"[1%]",tr_cvals[[3]],"[5%]",tr_cvals[[1]],"[10%]"), file="");

  # Max. eigenvalue test:
  vecm_eigen_coint <- ca.jo(zc, type="eigen", constant=FALSE, K=2, spec="longrun", season=NULL, dumvar=NULL, ctable="A2");
  max_statistics <- vecm_eigen_coint@teststat;
  max_cvals <- vecm_eigen_coint@cval;

  write(paste("\n* Max. eigenvalue statistic[e2] (H0: r=0 cointegrating vectors vs. H1: r=1) =", round(max_statistics[[2]],4)), file="");
  write(paste("* Max. eigenvalue statistic[e1] (H0: r=1 cointegrating vector vs. H1: r=2) =", round(max_statistics[[1]],4)), file="");
  write(paste("* Critical max-eigen test values [r=0]: ", max_cvals[[6]],"[1%]",max_cvals[[4]],"[5%]",max_cvals[[2]],"[10%]"), file="");
  write(paste("* Critical max-eigen test values [r=1]: ", max_cvals[[5]],"[1%]",max_cvals[[3]],"[5%]",max_cvals[[1]],"[10%]"), file="");

  if( (abs(tr_statistics[[2]]) < abs(tr_cvals[[1]])) && (abs(max_statistics[[2]]) < abs(max_cvals[[1]])) ) {
    write(paste("\n*** Series pair (xc,yc) is not cointegrating [H0: 0 cointegrating vectors accepted at >10% level]"),file="");
  } else if( (abs(tr_statistics[[2]]) < abs(tr_cvals[[3]])) && (abs(max_statistics[[2]]) < abs(max_cvals[[3]])) ) {
    write(paste("\n*** Series pair (xc,yc) is not cointegrating [H0: 0 cointegrating vectors accepted at >5% level]"),file="");
  } else if( (abs(tr_statistics[[1]]) < abs(tr_cvals[[1]])) && (abs(max_statistics[[1]]) < abs(max_cvals[[1]])) ) {
    write(paste("\n*** Series pair (xc,yc) is cointegrating [H0: 1 cointegrating vector accepted at >10% level]"),file="");
  } else if( (abs(tr_statistics[[1]]) < abs(tr_cvals[[3]])) && (abs(max_statistics[[1]]) < abs(max_cvals[[3]])) ) {
    write(paste("\n*** Series pair (xc,yc) is cointegrating [H0: 1 cointegrating vector accepted at >5% level]"),file="");
  } else if( (abs(tr_statistics[[1]]) > abs(tr_cvals[[5]])) && (abs(max_statistics[[1]]) > abs(max_cvals[[5]])) ) {
    write(paste("\n*** Series pair (xc,yc) is not cointegrating [H0: 1 cointegrating vector rejected at <1% level]"),file="");
  } else if( (abs(tr_statistics[[1]]) > abs(tr_cvals[[3]])) && (abs(max_statistics[[1]]) > abs(max_cvals[[3]])) ) {
    write(paste("\n*** Series pair (xc,yc) is not cointegrating [H0: 1 cointegrating vector rejected at <5% level]"),file="");
  }

} else {
  write(paste("\n*** Series pair (xc,yc) is not cointegrating since series are not integrated with same order."),file="");
}

#---- xnc, ync (non-cointegrated) series pair: ----

if( integ_flag_xnc && integ_flag_ync ) {
  # Trace test:
  vecm_trace_coint <- ca.jo(znc, type="trace", constant=FALSE, K=2, spec="longrun", season=NULL, dumvar=NULL, ctable="A2");
  eigenvalues <- vecm_trace_coint@lambda;
  tr_statistics <- vecm_trace_coint@teststat;
  tr_cvals <- vecm_trace_coint@cval;

  write(paste("* With max. lag = 2 and linear trend included."), file="");
  write(paste("* Eigenvalue (e1) =", round(eigenvalues[[2]],8),";","Eigenvalue (e2) =",round(eigenvalues[[1]],8)), file="");
  write(paste("\n* Trace statistic[e1 + e2] (H0: r=0 cointegrating vectors vs. H1: r=2) =", round(tr_statistics[[2]],4)), file="");
  write(paste("* Trace statistic[e1] (H0: r<=1 cointegrating vector vs. H1: r=2) =", round(tr_statistics[[1]],4)), file="");
  write(paste("* Critical trace test values [r=0]: ", tr_cvals[[6]],"[1%]",tr_cvals[[4]],"[5%]",tr_cvals[[2]],"[10%]"), file="");
  write(paste("* Critical trace test values [r<=1]: ", tr_cvals[[5]],"[1%]",tr_cvals[[3]],"[5%]",tr_cvals[[1]],"[10%]"), file="");

  # Max. eigenvalue test:
  vecm_eigen_coint <- ca.jo(znc, type="eigen", constant=FALSE, K=2, spec="longrun", season=NULL, dumvar=NULL, ctable="A2");
  max_statistics <- vecm_eigen_coint@teststat;
  max_cvals <- vecm_eigen_coint@cval;

  write(paste("\n* Max. eigenvalue statistic[e2] (H0: r=0 cointegrating vectors vs. H1: r=1) =", round(max_statistics[[2]],4)), file="");
  write(paste("* Max. eigenvalue statistic[e1] (H0: r=1 cointegrating vector vs. H1: r=2) =", round(max_statistics[[1]],4)), file="");
  write(paste("* Critical max-eigen test values [r=0]: ", max_cvals[[6]],"[1%]",max_cvals[[4]],"[5%]",max_cvals[[2]],"[10%]"), file="");
  write(paste("* Critical max-eigen test values [r=1]: ", max_cvals[[5]],"[1%]",max_cvals[[3]],"[5%]",max_cvals[[1]],"[10%]"), file="");

  if( (abs(tr_statistics[[2]]) < abs(tr_cvals[[1]])) && (abs(max_statistics[[2]]) < abs(max_cvals[[1]])) ) {
    write(paste("\n*** Series pair (xnc,ync) is not cointegrating [H0: 0 cointegrating vectors accepted at >10% level]"),file="");
  } else if( (abs(tr_statistics[[2]]) < abs(tr_cvals[[3]])) && (abs(max_statistics[[2]]) < abs(max_cvals[[3]])) ) {
    write(paste("\n*** Series pair (xnc,ync) is not cointegrating [H0: 0 cointegrating vectors accepted at >5% level]"),file="");
  } else if( (abs(tr_statistics[[1]]) < abs(tr_cvals[[1]])) && (abs(max_statistics[[1]]) < abs(max_cvals[[1]])) ) {
    write(paste("\n*** Series pair (xnc,ync) is cointegrating [H0: 1 cointegrating vector accepted at >10% level]"),file="");
  } else if( (abs(tr_statistics[[1]]) < abs(tr_cvals[[3]])) && (abs(max_statistics[[1]]) < abs(max_cvals[[3]])) ) {
    write(paste("\n*** Series pair (xnc,ync) is cointegrating [H0: 1 cointegrating vector accepted at >5% level]"),file="");
  } else if( (abs(tr_statistics[[1]]) > abs(tr_cvals[[5]])) && (abs(max_statistics[[1]]) > abs(max_cvals[[5]])) ) {
    write(paste("\n*** Series pair (xnc,ync) is not cointegrating [H0: 1 cointegrating vector rejected at <1% level]"),file="");
  } else if( (abs(tr_statistics[[1]]) > abs(tr_cvals[[3]])) && (abs(max_statistics[[1]]) > abs(max_cvals[[3]])) ) {
    write(paste("\n*** Series pair (xnc,ync) is not cointegrating [H0: 1 cointegrating vector rejected at <5% level]"),file="");
  }

} else {
  write(paste("*** Series pair (xnc,ync) is not cointegrating since series are not integrated with same order."),file="");
}

dev.off()
q()