TEXT   99

# 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.
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. }