cc <- function(y, adf=TRUE, dtt=TRUE, graph=TRUE, pvar=4, ira=24){
  if(!interactive())
    {
      return("The function can only be executed in interactive mode.")
    }
  Tall <- nrow(y)
  Tsub <- Tall-pvar
  K <- ncol(y)
  pcc <- pvar-1
  ly <- rbind(rep(NA,1),y)[1:Tall,]
  dy <- y-ly
  if(Tsub<K*(pcc+1))
  {
    return("Insufficient degrees of freedom.")
  }
  if(pvar<=0)
  {
    return("Please provide a positive integer value for the lag length of the underlying VAR-process.")
  }
#
  if(adf)
  {
    library(tseries)
    print(apply(y,2,adf.test))
    print(apply(na.omit(dy),2,adf.test))
  }
#  
  if(pcc >= 1)
  {
    dyl <- array(NA,c(Tall,K,pcc))
    dat <- dy
    for(i in 1:pcc)
    {
      nval <- matrix(NA,i,K)
      dyl[,,i] <- rbind(nval,dy)[1:Tall,]
      dat <- cbind(dat,dyl[,,i])
    }
    dattsub <- dat[-(1:pvar),-(1:K)]
    dytsub <- dy[-(1:pvar),]
    lytsub <- ly[-(1:pvar),]
    if(dtt)
    {
      uhat <- lm(dytsub~dattsub)$resid
      ucoef <- lm(dytsub~dattsub)$coefficients
      vhat <- lm(lytsub~dattsub)$resid
      vcoef <- lm(lytsub~dattsub)$coefficients
    }
    else
    {
      one <- matrix(1,Tsub,1)
      uhat <- lm(dytsub~-1+dattsub)$resid
      ucoef <- lm(dytsub~-1+dattsub)$coefficients
      w1hat <- lm(one~-1+dattsub)$resid
      w1coef <- lm(one~-1+dattsub)$coefficients
      w2hat <- lm(lytsub~-1+dattsub)$resid
      w2coef <- lm(lytsub~-1+dattsub)$coefficients
      vhat <- cbind(w1hat,w2hat)
    }
    sig.uhat <- 1/Tsub*(t(uhat)%*%uhat)
    sig.vhat <- 1/Tsub*(t(vhat)%*%vhat)
    sig.uvhat <- 1/Tsub*(t(uhat)%*%vhat)
  }
  else
  {
    dytsub <- dy[-1,]
    lytsub <- ly[-1,]
    if(dtt)
      {
       uhat <- dytsub-apply(dytsub,2,mean)
       vhat <- lytsub-apply(lytsub,2,mean)
      }
    else
      {
        one <- matrix(1,Tsub,1)
        uhat <- dytsub
        w1hat <- one
        w2hat <- lytsub
        vhat <- cbind(w1hat,w2hat)
      }
    sig.uhat <- 1/Tsub*(t(uhat)%*%uhat)
    sig.vhat <- 1/Tsub*(t(vhat)%*%vhat)
    sig.uvhat <- 1/Tsub*(t(uhat)%*%vhat)
  }
  M <- solve(sig.vhat) %*% t(sig.uvhat) %*% solve(sig.uhat) %*% sig.uvhat
  ev <- eigen(M)$values
  evlr <- -Tsub*log(1-ev)
  print(evlr)
  if(dtt)
    {
      trans <- upper.tri(matrix(1,K,K),diag=T)
      tr <- trans%*%evlr
    }
  else
    {
     trans <- upper.tri(matrix(1,K+1,K+1),diag=T)
     tr <- trans%*%evlr
    } 
  print(tr)
  print("Please, provide the number of cointegrating vectors.")
  h <- as.integer(abs(scan(file="", nmax=1)))
  if(h<1 || h>length(tr))
   {
    print("Please provide a positive integer value that is neither smaller than 1 nor greater than the maximum number of cointegrating vectors.")
    h <- as.integer(abs(scan(file="", nmax=1)))
   }
  Aorig <- eigen(M)$vectors
  A2 <- colSums(Aorig*Aorig)
  Anorm <- sweep(Aorig,2,STATS=sqrt(A2),FUN="/")
  norm <- function(X) X / sqrt(t(X)%*%sig.vhat%*%X)
  Atmp <- apply(Anorm,2, norm)
  A <- matrix(Atmp,length(A2),length(A2))
  A <- A[,1:h]
  if(dtt)
    {
      if(pcc==0)
        {
          Psi0 <<- sig.uvhat%*%A%*%t(A)
          alpha <<- apply(dytsub,2,mean)-Psi0%*%apply(lytsub,2,mean)
          Omega <- 1/Tsub*t(uhat-vhat%*%Psi0)%*%(uhat-vhat%*%Psi0)
          Psi <<- "Not applicable - no lags specified."
          print(alpha)
          print(Psi0)
          print(Omega)
        }
      else
        {
          Psi0 <- sig.uvhat%*%A%*%t(A)
          Omega <- 1/Tsub*t(uhat-vhat%*%Psi0)%*%(uhat-vhat%*%Psi0)
          pi0 <- ucoef[1,]
          xi0 <- vcoef[1,]
          alpha <- pi0-Psi0%*%xi0
          Piall <- ucoef[-1,]
          Xiall <- vcoef[-1,]
          ind1 <- seq(1,K*pcc,K)
          ind2 <- seq(K,K*pcc,K)
          Pi <- array(NA,c(K,K,pcc))
          Xi <- array(NA,c(K,K,pcc))
          Psi <- array(NA,c(K,K,pcc))
          for(i in 1:pcc)
            {
              Pi[,,i] <- Piall[ind1[i]:ind2[i],]
              Xi[,,i] <- Xiall[ind1[i]:ind2[i],]
              Psi[,,i] <- Pi[,,i]- Psi0%*%Xi[,,i]
            }
          print(alpha)
          print(Psi0)
          print(Psi)
          print(Omega)
        }
    }
  else
    {
      if(pcc==0)
        {
          Psi0 <- sig.uvhat%*%A%*%t(A)[,-1]
          alpha <- sig.uvhat%*%A%*%t(A)[,1]
          Omega <- 1/Tsub*t(uhat-w1hat%*%t(alpha)-w2hat%*%Psi0)%*%(uhat-w1hat%*%t(alpha)-w2hat%*%Psi0)
          Psi <- "Not applicable - no lags specified."
          print(alpha)
          print(Psi0)
          print(Omega)
        }
      else
        {
          Psi0 <- sig.uvhat%*%A%*%t(A)[,-1]
          alpha <- sig.uvhat%*%A%*%t(A)[,1]
          Omega <- 1/Tsub*t(uhat-w1hat%*%t(alpha)-w2hat%*%Psi0)%*%(uhat-w1hat%*%t(alpha)-w2hat%*%Psi0)
          Piall <- ucoef
          Nyall <- w1coef
          Xiall <- w2coef
          ind1 <- seq(1,K*pcc,K)
          ind2 <- seq(K,K*pcc,K)
          Pi <- array(NA,c(K,K,pcc))
          Ny <- array(NA,c(K,1,pcc))
          Xi <- array(NA,c(K,K,pcc))
          Psi <- array(NA,c(K,K,pcc))
          for(i in 1:pcc)
            {
              Pi[,,i] <- Piall[ind1[i]:ind2[i],]
              Ny[,1,i] <- Nyall[ind1[i]:ind2[i]]
              Xi[,,i] <- Xiall[ind1[i]:ind2[i],]
              Psi[,,i] <- Pi[,,i]-alpha%*%t(Ny[,,i])-Psi0%*%Xi[,,i]
            }
          print(alpha)
          print(Psi0)
          print(Psi)
          print(Omega)
        }
    }
  endog <- dytsub
  datmat <- matrix(NA,K,K)
  const <- matrix(1,Tsub,1)
  fit <- matrix(NA,Tsub,K)
  if(pcc==0)
   {
    fittedd <- const%*%t(alpha)+lytsub%*%t(Psi0)
    varc <- ar.ols(y, aic=F, order.max=1)$ar
    Psi0L <- varc[1,,]
    fittedl <- const%*%t(alpha)+lytsub%*%t(Psi0L)
   }
  if(pcc>0)
   {
    datmat <- Psi[,,1]
   }
  if(pcc>1)
   {
    for(i in 2:pcc)
     {
      datmat <- cbind(datmat,Psi[,,i])
     }
  fittedd <- const%*%t(alpha)+dattsub%*%t(datmat)+lytsub%*%t(Psi0)
  psum <- function(X)
    {
     tmp <- matrix(0,K,K)
     for(i in 1:pvar)
      {
       tmp <- tmp+X[i,,]
      }
    }
  varc <- ar.ols(y, aic=F, order.max=pvar)$ar
  Psi0L <- psum(varc)
  fittedl <- const%*%t(alpha)+dattsub%*%t(datmat)+lytsub%*%t(Psi0L)
  }
  if(graph)
    if (file.access("graph_y1_L.jpg")==0)
      {
        system(paste("cmd /c cd ", getwd(), " && del /q *.jpg", sep=""))
      }
  {
    for (i in 1:K)
      {
        jpeg(filename=paste("graph_y",i,"_L.jpg",sep=""), width = 480, height = 480)
        ts.plot(ts(y[-(1:pvar),i]),ts(fittedl[,i]), col=c("black","red"),main=substring(paste('Diagram of fit: y',i,' (levels)',sep=""),1,60),xlab="Index")
        dev.off()
        jpeg(filename=paste("graph_y",i,"_D.jpg",sep=""), width = 480, height = 480)
        ts.plot(ts(dytsub[,i]),ts(fittedd[,i]), col=c("black","red"),main=substring(paste('Diagram of fit: y',i,' (1st differences)',sep=""),1,60),xlab="Index")
        dev.off()
      }
  }
# Impulse Response Functions
  ID <- diag(K)
  NM <- matrix(0,nrow=K,ncol=K*pvar)
  D <- array(NA, c(K,K,pvar))
  B <- array(NA, c(K,K,ira))
  dysub <- dy[-(1:pvar),]
  RES <- dysub-fittedd
  SIG <- var(RES)
  GAM <- t(chol(SIG))
  GAM <- lower.tri(GAM)*GAM+ID
  if(pvar==1)
    {
      D[,,1] <- ID-Psi0
    }
  if(pvar==2)
    {
      D[,,1] <- ID-Psi0+Psi[,,1]
      D[,,2] <- -Psi[,,pcc]
    }
  if(pvar>2)
    {
      D[,,1] <- ID-Psi0+Psi[,,1]
      for (i in 2:pcc)
        {
          D[,,i] <- Psi[,,i]-Psi[,,i-1]
        }
      D[,,pvar] <- -Psi[,,pcc]
    }
  Dstack <- D[,,1]
  if(pvar>1)
    {
      for(i in 2:pvar)
        {
          Dstack <- rbind(Dstack,D[,,i])
        }
    }
  tmp <- cbind(ID,NM)[1:K,1:(K*pvar)]
  for(i in 1:ira)
    {
      tmpB <- tmp%*%Dstack
      tmp <- cbind(tmpB,tmp)[1:K,1:(K*pvar)]
      B[,,i] <- tmpB%*%GAM
    }
  if(graph)
  {
    for(i in 1:K)
      {
        for(j in 1:K)
          {
            jpeg(filename=paste("ira_y",i,"_on_y",j,".jpg",sep=""))
           barplot(B[j,i,],col="seagreen",names.arg=c(1:ira),main="Cumulated, Orthogonal Impulse Response",sub=paste("Unit shock in y",i," on y",j, sep="" ))
           dev.off()
         }
      }
  }
# End of Impulse Response Function
  structure(list(ALPHA=alpha, PSI0=Psi0,PSI=Psi,OMEGA=Omega, ENDOG=endog, FITTEDD=fittedd, FITTEDL=fittedl, IRA=B))
}
