## (a) sample variances, along the main diagonal,
## and intervenor-adjusted partial correlations, below the main diagonal

Stats.ia<-function(data){
  n <- ncol(data)
  V <- matrix(rep(0,n*n), nrow = n)
  var <- as.vector(diag(cov(data)))
  cor <- as.matrix(cor(data))     
  for(i in 1:(n-1)){

    ## sample variances on main diagonal 
    V[i,i] <- var[i]

    ## sample correlations on the first off-diagonal              
    V[(i+1),i] <- cor[(i+1),i]    
  }
  V[n,n] <- var[n]      
  for(i in 3:n){
    for(j in 1:(i-2)){
      A <- cbind(data[,i],data[,j])
      B <- as.matrix(data[,(j+1):(i-1)])

      ## sample intervenor-adjusted partial correlations on the pth off-diagonal(p=2...n-1)
      V[i,j] <- cor(residuals(lm(A~B)))[1,2]  
    }
  } 
  V
}


## (b) sample partial variances, along the main diagonal,
## and partial correlations, below the main diagonal

Stats.pc<-function(data){
  n <- ncol(data)
  V <- matrix(rep(0,n*n), nrow = n)
  for(i in 1:n){
    for(j in 1:i){
      if(i==j){
        A <- data[,i]
        B <- as.matrix(data[,-i])

        ## sample partial variances on the main diagonal  
        V[i,i] <- var(residuals(lm(A~B)))      
      }            
      else{   
        A <- cbind(data[,i],data[,j])
        B <- as.matrix(data[,-c(i,j)])

        ## sample partial correlations below the main diagonal                
        V[i,j] <- cor(residuals(lm(A~B)))[1,2] 
      }
    }  
  } 
  V
}


## (c) sample innovation variances, along the main diagonal,
## and autoregressive coefficients, below the main diagonal

Stats.ar<-function(data, MSE){
  n <- ncol(data)
  m <- nrow(data)
  V <- matrix(rep(0,n*n), nrow = n)
  cov <- matrix(cov(data), nrow=n) 
  V[1,1] <- cov[1,1]
  for(i in 2:n){
    A <- data[,i]
    B <- as.matrix(data[,1:(i-1)])

    ## innovation variances_residual variance(2.18)
    if(MSE==TRUE)  V[i,i] <- var(residuals(lm(A~B)))*(m-1)/(m-i)  
    else  V[i,i] <- var(residuals(lm(A~B)))
      
    ## autoregressive coefficients(2.19)         
    phi <- as.vector(t(cov[1:(i-1),i])%*%solve(cov[1:(i-1),1:(i-1)]))   
    for(j in 1:(i-1)){
      V[i,j] <- phi[j]
    }
  } 
  V
}


## Conclusion for Stats.ia
flag.ia<-function(data){
  N <- nrow(data)
  n <- ncol(data)
  r.thumb <- matrix(rep(0, n*n), nrow = n)
  V <- Stats.ia(data)

  for(i in 2:n){
    for(j in 1:(i-1)){

      c <- 2/sqrt(N-(i-j-1))
      if(abs(V[i,j]) > c)   r.thumb[i, j] <- "r"
      else   r.thumb[i, j] <- "a"
    }
  } 
  list(flag.ia = r.thumb)    
}


## Conclusion for Stats.pc
flag.pc<-function(data){
  N <- nrow(data)
  n <- ncol(data)
  r.thumb <- matrix(rep(0, n*n), nrow = n)

  V <- Stats.pc(data)
  for(i in 2:n){
    for(j in 1:(i-1)){

        c <- 2/sqrt(N-(n-2))
      if(abs(V[i,j]) > c)   r.thumb[i, j] <- "r"
      else   r.thumb[i, j] <- "a"
    }
  } 
  list(flag.pc = r.thumb)    
}

## Conclusion for Stats.ar (the code is the same as the code for theorem6.3(a))