R/disag1.R

In disagg2: Support Functions for Time Series Analysis Book

mySym <- function(x) {
n1 <- length(x)
y <- matrix(0,nrow=n1,ncol=n1)
j <- n1
for(i in 1:n1) {
   y[i,i:n1] <- x[1:j]
   j <- j - 1
   }
   y[lower.tri(y)] <- t(y)[lower.tri(y)]
   return(y)
   }


disag1 <- function(x,m) {
if(tsp(x)[3]*m!=12)stop("freq and m don't work")
y1 <- mean(x)
x1 <- x - y1
y <- arima(x1,order=c(1,0,1),include.mean=FALSE)
N <- length(x)
n1 <- N * m
n3 <- (N*m) + (m-1)

gammau <- length(N)
gammaw <- length(n1)
gammau[1] <- ((1 - 2*y$coef[1]*y$coef[2] + (y$coef[2]^2)) * y$sigma2)/(1 - y$coef[1]^2)
gammau[2] <- y$coef[1]*gammau[1] - (y$coef[2]*y$sigma2)
for(i in 3:N) {
gammau[i] <- y$coef[1]*gammau[(i-1)]
}
p1 <- polynom(rep(1,m))^2
c1 <- coef(p1)
n2 <- length(c1)
p3 <- polynom(rep(1,m))
c3 <- coef(p3)
nc3 <- length(c3)
c2 <- c(m,2*((m-1):1))
amat <- matrix(0,nrow=2,ncol=(2*m))
amat[1,1:m] <- c2
amat[2,2:(2*m)] <- c1
gmat <- matrix(0,nrow=(2*m),ncol=2)
gmat[1,1] <- 1
gmat[2,2] <- 1

cmat <- matrix(0,nrow=N,ncol=n1)
for(i in 1:N) {
j <- (i-1)*m + 1
cmat[i,j:(j+nc3-1)]<- c3
}
xs <- sign(y$coef[1])
xp <- abs(y$coef[1])^(1/m)
newphi <- xp * xs
phiv <- newphi^(1:(m+1))
gmat[m:(2*m),2] <- phiv
biga <- amat %*% gmat
gs <- solve(biga,matrix(gammau[1:2],ncol=1))
gammaw[1:2] <- gs
for(i in 3:n1) {
gammaw[i] <- newphi*gammaw[(i-1)]
}
gwmat <- mySym(gammaw)
gumat <- mySym(gammau)
gumati <- solve(gumat)
gy <- gumati %*% x
bigy_s <- gwmat %*% t(cmat) %*% gy
bigy_m <- bigy_s * m
bigy_ar <- arima(bigy_s,order=c(1,0,1))
zzz <- list(y_s=bigy_s,y_m=bigy_m,disphi=newphi,distheta=bigy_ar$coef[2],
    dissig2=bigy_ar$sigma2)
return(zzz)
}