TEXT   83

The following functions only work for constant order AD models

Guest on 23rd April 2022 01:13:09 AM

  1.  
  2. ##The following functions only work for constant order AD models
  3. ###############Theorem 5.8 (a)###############  
  4.  
  5. REMLE.phi.delta<-function(data, covariate, p, REMLE){
  6.   data <- as.matrix(data)
  7.   Z <- as.matrix(covariate)
  8.   n <- ncol(data)
  9.   m <- ncol(Z)
  10.   V <- matrix(rep(0,n*n), nrow = n)
  11.  
  12.   if(length(p)==1)  p <- rep(0:p, c(rep(1,p),(n-p)))  
  13.  
  14.   for(i in 1:n){
  15.     Z <- as.matrix(Z[!is.na(data[,i]),])    
  16.     data <- data[!is.na(data[,i]),]          
  17.     N <- nrow(data)        
  18.  
  19.  
  20.     if(m==1)  A <- cov(matrix(data[ ,(1:i)], ncol=i))*(N-m)/N  
  21.          
  22.     else{
  23.       Y <- NULL
  24.       for(j in 1:N)
  25.         Y <- cbind(Y, t(data[j,1:i]))
  26.       Y <- t(Y)
  27.  
  28.       ## p.147
  29.       beta.bar <- kronecker(solve(t(Z)%*%Z)%*%t(Z), diag(i))%*%Y
  30.       A <- matrix(rep(0,i*i), nrow = i)
  31.       for(s in 1:N){
  32.         M <- matrix(data[s,1:i], nr=i)-kronecker(matrix(Z[s,], nr=1), diag(i))%*%beta.bar
  33.         A <- A + M %*% t(M)
  34.       }
  35.       A <- A/N
  36.     }
  37.  
  38.     if(REMLE==TRUE)  cov <- (N/(N-m))*A      
  39.     else  cov <- A
  40.  
  41.  
  42.     ## innovation variances(2.18) for i=1
  43.     if(i==1)  V[i,i] <- cov
  44.  
  45.     else{
  46.  
  47.       ## innovation variances(2.18)
  48.       V[i,i] <- cov[i,i]-t(cov[(i-min(p[i],(i-1))):(i-1), i]) %*%
  49.          solve(cov[(i-min(p[i],(i-1))):(i-1), (i-min(p[i],(i-1))):(i-1)]) %*%
  50.          cov[(i-min(p[i],(i-1))):(i-1), i]    
  51.  
  52.       ## autoregressive coefficients(2.19)
  53.       phi <- as.vector(t(cov[(i-min(p[i],(i-1))):(i-1),i]) %*%
  54.          solve(cov[(i-min(p[i],(i-1))):(i-1),(i-min(p[i],(i-1))):(i-1)]))  
  55.       for(j in (i-min(p[i],(i-1))):(i-1))
  56.         V[i,j] <- phi[j+1-(i-min(p[i],(i-1)))]
  57.     }
  58.   }  
  59.   V
  60. }
  61.  
  62.  
  63. ###############Theorem 5.8 (b)###############  
  64.  
  65. ## calculate the correlations below and above the main diagonals
  66. cov.cor<-function(V){
  67.   for(i in 2: nrow(V)){
  68.     for(j in 1:(i-1))
  69.       V[j,i] <- V[i,j] <- V[i,j]/sqrt(V[i,i]*V[j,j])
  70.   }
  71.   V
  72. }
  73.  
  74.  
  75. REMLE.cov<-function(data, covariate, p, REMLE){
  76.   n <- ncol(data)
  77.   T <- -REMLE.phi.delta(data, covariate, p, REMLE=REMLE)
  78.   D <- matrix(rep(0,n*n), nrow = n)
  79.   diag(D) <- -diag(T)
  80.   diag(T) <- rep(1,n)
  81.  
  82.   solve(T) %*% D %*% solve(t(T))
  83. }
  84.  
  85.  
  86. ###############Theorem 5.8 (c)###############  
  87. ## inverse
  88. inv <- function(cov){
  89.   solve(cov)
  90. }
  91.  
  92.  
  93. ## intervenor-adjusted partial correlation
  94.  
  95. ## recursion formula for partial correlations in (2.6)
  96. recur.pcor<-function(index1, index2, index3, V){
  97.   index1 <- index1
  98.   index2 <- index2
  99.   index3 <- index3
  100.   index <- c(index1, index2, index3)
  101.  
  102.   if(length(index3)==1){
  103.     rxy.z <- (V[index[1], index[2]]-V[index[1], index[3]]*V[index[2], index[3]])/
  104.             (sqrt(1-V[index[1], index[3]]^2)*sqrt(1-V[index[2], index[3]]^2))
  105.     return(rxy.z)
  106.   }
  107.   else{
  108.     index1 <- index[1]
  109.     index2 <- index[2]
  110.     index30 <- index[3]    
  111.     index3c <- index[-c(1,2,3)]
  112.  
  113.     rxy.zc <- recur.pcor(index1, index2, index3c, V)
  114.     rxz0.zc <- recur.pcor(index1,index30,index3c, V)
  115.     ryz0.zc <- recur.pcor(index2,index30,index3c, V)
  116.     rxy.z <- (rxy.zc - rxz0.zc*ryz0.zc)/(sqrt(1-rxz0.zc^2)*sqrt(1-ryz0.zc^2))
  117.     return(rxy.z)
  118.   }                      
  119. }  
  120.  
  121. partial.inter<-function(cor){
  122.   V <- cor
  123.   n <- ncol(cor)
  124.   m <- n-2
  125.   for(k in 1:m){
  126.     if(k==m)
  127.        V[1, n] <- V[n, 1] <- recur.pcor(1, n, c(2:(n-1)), cor)
  128.  
  129.      else
  130.        diag(V[1:(n-k-1), (k+2):n]) <- diag(V[(k+2):n, 1:(n-k-1)]) <- sapply(1:(m+1-k), function(i) if(k==1)
  131.                     recur.pcor(i, i+k+1, i+1, cor) else recur.pcor(i, i+k+1, c((i+1):(i+k)), cor))
  132.   }
  133.   V
  134. }
  135.  
  136.  
  137.  
  138. ###############Theorem 5.8 (d)###############
  139.  
  140. beta.hat<-function(data, covariate, p){
  141.   data <- as.matrix(data)
  142.   Z <- as.matrix(covariate)
  143.   n <- ncol(data)
  144.   m <- ncol(covariate)
  145.  
  146.   if(length(p)==1)  p <- rep(0:p, c(rep(1,p),(n-p)))  
  147.  
  148.   T <- -REMLE.phi.delta(data, covariate, p, REMLE=TRUE)  ## REMLE=TRUE or FALSE
  149.   diag(T) <- rep(1,n)
  150.  
  151.   if(m==1){
  152.  
  153.     muhatstar <- rep(0,n)
  154.     for(i in 1:n){  
  155.       data1 <- data[!is.na(data[,i]),]      
  156.  
  157.       if(i==1)  muhatstar[i] <- mean(data1[ ,1])
  158.       else{
  159.         Y <- matrix(apply(data1[ ,(1:i)], 2, mean),nrow=1)
  160.         muhatstar[i] <- Y[i]+T[i,(i-min(p[i],(i-1))):(i-1)] %*% Y[(i-min(p[i],(i-1))):(i-1)]
  161.       }
  162.     }
  163.     muhat <- solve(T) %*% muhatstar
  164.     list(muhat=muhat)
  165.   }
  166.  
  167.   else{
  168.     betahatstar <- NULL
  169.     for(i in 1:n){  
  170.       Z1 <- Z[!is.na(data[,i]),]
  171.       data1 <- data[!is.na(data[,i]),]
  172.       N1 <- nrow(data1)
  173.  
  174.       if(min(p[i],(i-1))==0)  betahatstar <- rbind(betahatstar, solve(t(Z1) %*% Z1) %*% t(Z1) %*% data1[, i])
  175.       else{
  176.         Q <- NULL
  177.         for(k in 1:min(p[i],(i-1)))
  178.           Q <- cbind(Q, data1[, i-k])
  179.         Q <- matrix(Q, nrow=N1)
  180.  
  181.         M <- diag(N1) - Q %*% solve(t(Q) %*% Q) %*% t(Q)  
  182.         betahatstar <- rbind(betahatstar, solve(t(Z1) %*% M %*% Z1) %*% t(Z1) %*% M %*% data1[, i])
  183.       }
  184.     }
  185.     betahat <- kronecker(solve(T),diag(m)) %*% betahatstar
  186.     list(betahat=betahat)
  187.   }
  188. }

Raw Paste

Login or Register to edit or fork this paste. It's free.