# R/mvrnormSeries.R In sae2: Small Area Estimation: Time-Series Models

#### Documented in mvrnormSeries

```mvrnormSeries <- function(NV=1, D, T, sigma.e,
rho.dyn, sigma.v.dyn, sigma.u.dyn, rho.u.dyn,
rho.RY, sigma.v.RY, sigma.u.RY, rho.u.RY,
tol=1e-6, conditional.mean=FALSE) {
if (missing(T)) stop("'T' must be specified")
if (T < 3) stop("'T' must be >2")
if (class(sigma.e) == "list") {
if (missing(D)) {
D <- length(sigma.e)
} else {
if (D != length(sigma.e))
stop("the length of 'sigma.e' must agree with D")
}
lapply(sigma.e, FUN = function(x) { if (!all(dim(x) == c(NV*T, NV*T) ))
stop(paste("each element of 'sigma.e' must be a square matrix",
"with 'NV*T' rows"))
} )
} else {
if (missing(D)) stop("'D' must be specified")
if (class(sigma.e) != "matrix")
stop("'sigma.e' must be a matrix or a list of matrices")
if (!all(dim(sigma.e) == c(NV*T, NV*T) ))
stop("'sigma.e' must be a square matrix with 'NV*T' rows")
}
sigma <- matrix(0, NV*T, NV*T)
mu <- rep(0, NV*T)
if (!(missing(rho.dyn) | missing(sigma.v.dyn) | missing(sigma.u.dyn))) {
if (missing(rho.dyn) | missing(sigma.v.dyn) | missing(sigma.u.dyn))
stop("some parameters missing for dynamic model")
if (NV ==1) {
if(length(sigma.v.dyn) != 1 |
length(sigma.u.dyn) != 1)
stop("wrong length for 'sigma.v.dyn' or 'sigma.u.dyn'")
sigma <- sigma.v.dyn*Gamma_v_f(rho.dyn, T)[[1]] +
sigma.u.dyn*Gamma_u_f(rho.dyn, T)[[1]]
} else {
if(missing(rho.u.dyn)) stop("'rho.u.dyn' must be specified")
len.rho <- (NV*(NV-1))/2
if(length(rho.u.dyn)==1 & NV>2) rho.u.dyn <- rep(rho.u.dyn, len.rho)
if(length(rho.u.dyn)!= len.rho)
stop(paste0("the length of 'rho.u.dyn' must be 1 or (NV*(NV+1))/2"))
u_corr <- diag(1, NV)
for (i in 2:NV) for (j in 1:(i-1)){
u_corr[i,j] <- u_corr[j,i] <- rho.u.dyn[((i-2)*(i-1))/2+j]
}
if(length(sigma.v.dyn)==1) sigma.v.dyn <- rep(sigma.v.dyn, NV)
if(length(sigma.u.dyn)==1) sigma.u.dyn <- rep(sigma.u.dyn, NV)
if(length(sigma.v.dyn)!= NV)
stop(paste0("the length of 'sigma.v.dyn' must be 1 or NV"))
if(length(sigma.u.dyn)!= NV)
stop(paste0("the length of 'sigma.u.dyn' must be 1 or NV"))
sqrt_u <- sqrt(sigma.u.dyn)
sqrt_v <- sqrt(sigma.v.dyn)
sigma  <-
((sqrt_u%*%t(sqrt_u)) * u_corr) %x% Gamma_u_f(rho.dyn, T)[[1]] +
((sqrt_v%*%t(sqrt_v)) * u_corr) %x% Gamma_v_f(rho.dyn, T)[[1]]
}
}
if (!(missing(rho.RY) | missing(sigma.v.RY) | missing(sigma.u.RY))) {
if (missing(rho.RY) | missing(sigma.v.RY) | missing(sigma.u.RY))
stop("some parameters missing for Rao-Yu model")
Gamma_v <- matrix(1,nrow=T,ncol=T)
Gamma_b <- matrix(rep(1:T,times=T),nrow=T, ncol=T)
Gamma_s <- abs(Gamma_b-t(Gamma_b))
if(rho.RY > 0) {
Gamma_u <- rho.RY^(Gamma_s)/(1-rho.RY^2)
} else {
Gamma_u <- diag(T)
}
if (NV ==1) {
if(length(sigma.v.RY) != 1 |
length(sigma.u.RY) != 1)
stop("wrong length for 'sigma.v.RY' or 'sigma.u.RY'")
sigma <- sigma + sigma.v.RY*Gamma_v +
sigma.u.RY*Gamma_u
} else {
if(missing(rho.u.RY)) stop("'rho.u.RY' must be specified")
len.rho <- (NV*(NV-1))/2
if(length(rho.u.RY)==1 & NV>2) rho.u.RY <- rep(rho.u.RY, len.rho)
if(length(rho.u.RY)!= len.rho)
stop(paste0("the length of 'rho.u.RY' must be 1 or (NV*(NV+1))/2"))
u_corr <- diag(1, NV)
for (i in 2:NV) for (j in 1:(i-1)){
u_corr[i,j] <- u_corr[j,i] <- rho.u.RY[((i-2)*(i-1))/2+j]
}
if(length(sigma.v.RY)==1) sigma.v.RY <- rep(sigma.v.RY, NV)
if(length(sigma.u.RY)==1) sigma.u.RY <- rep(sigma.u.RY, NV)
if(length(sigma.v.RY)!= NV)
stop(paste0("the length of 'sigma.v.RY' must be 1 or NV"))
if(length(sigma.u.RY)!= NV)
stop(paste0("the length of 'sigma.u.RY' must be 1 or NV"))
sqrt_u <- sqrt(sigma.u.RY)
sqrt_v <- sqrt(sigma.v.RY)
sigma  <- sigma +
((sqrt_u%*%t(sqrt_u)) * u_corr) %x% Gamma_u +
((sqrt_v%*%t(sqrt_v)) * u_corr) %x% Gamma_v
}
}
if (NV == 1) {
y <- rep(0, D*T)
if(conditional.mean) {
y.mean <- rep(0, D*T)
for (i in 1:D) {
y.mean[((i-1)*T + 1):(i*T)] <-
mvrnorm(mu=mu, Sigma=sigma, tol=tol)
}
}
if(class(sigma.e) == "list") {
if(conditional.mean) {
for (i in 1:D) {
y[((i-1)*T + 1):(i*T)] <- y.mean[((i-1)*T + 1):(i*T)] +
mvrnorm(mu=mu, Sigma=sigma.e[[i]], tol=tol)
}
} else {
for (i in 1:D) {
y[((i-1)*T + 1):(i*T)] <-
mvrnorm(mu=mu, Sigma=sigma+sigma.e[[i]], tol=tol)
}
}
} else {
if(conditional.mean) {
for (i in 1:D) {
y[((i-1)*T + 1):(i*T)] <- y.mean[((i-1)*T + 1):(i*T)] +
mvrnorm(mu=mu, Sigma=sigma.e, tol=tol)
}
} else {
for (i in 1:D) {
y[((i-1)*T + 1):(i*T)] <-
mvrnorm(mu=mu, Sigma=sigma+sigma.e, tol=tol)
}
}
}
} else {
y <- matrix(0, nrow=D*T, ncol=NV)
if(conditional.mean) {
y.mean <- matrix(0, nrow=D*T, ncol=NV)
for (i in 1:D) {
y.mean.temp <- mvrnorm(mu=mu, Sigma=sigma, tol=tol)
for (nv in 1:NV) {
y.mean[((i-1)*T+1):(i*T), nv] <- y.mean.temp[((nv-1)*T+1):(nv*T)]
}
}
}
if(class(sigma.e) == "list") {
if(conditional.mean) {
for (i in 1:D) {
y.temp <- mvrnorm(mu=mu, Sigma=sigma.e[[i]], tol=tol)
for (nv in 1:NV) {
y[((i-1)*T+1):(i*T), nv] <- y.mean[((i-1)*T+1):(i*T), nv] +
y.temp[((nv-1)*T+1):(nv*T)]
}
}
} else {
for (i in 1:D) {
y.temp <-
mvrnorm(mu=mu, Sigma=sigma+sigma.e[[i]], tol=tol)
for (nv in 1:NV) {
y[((i-1)*T+1):(i*T), nv] <- y.temp[((nv-1)*T+1):(nv*T)]
}
}
}
} else {
if(conditional.mean) {
for (i in 1:D) {
y.temp <- mvrnorm(mu=mu, Sigma=sigma.e, tol=tol)
for (nv in 1:NV) {
y[((i-1)*T+1):(i*T), nv] <- y.mean[((i-1)*T+1):(i*T), nv] +
y.temp[((nv-1)*T+1):(nv*T)]
}
}
} else {
for (i in 1:D) {
y.temp <-
mvrnorm(mu=mu, Sigma=sigma+sigma.e, tol=tol)
for (nv in 1:NV) {
y[((i-1)*T+1):(i*T), nv] <- y.temp[((nv-1)*T+1):(nv*T)]
}
}
}
}
}
if(conditional.mean) {
return(list(y=y, y.mean=y.mean))
} else {
return(y)
}
}
```

## Try the sae2 package in your browser

Any scripts or data that you put into this service are public.

sae2 documentation built on July 22, 2019, 9:04 a.m.