Dear Splus users,
I do need some advices about vectors definition and storage.
I am building a quite big procedure where many vectors and matrices are
involved. Within this procedure there is a for loop (of length n),
and because of at each step I want to store all the vectors and matrices, I
use the following structure:
array(NA, dim = c(n, number of rows, 1)) or array(NA, dim = c(number of rows,
1)) for the vectors and
array(NA, dim = c(n, number of rows, number of columns)) or array(NA, dim =
c(number of rows, number of columns)) for the matrices.
For some reasons it happens that vectors are stored not correctly:
not vertically (for example of dimension 13 by 1)
but horizontally (i.e. of dimension 1 by 13).
In the function below reported, the vector E4iniz is correctly evaluated
(dimension (14 by 1), but with
E6iniz <- E4iniz[(2:14), ]
something goes wrong, E6iniz becomes a vector of 1 row and 13 columns.
Why this behaviour? Have you got any hint or suggestion?
The input parameters for the procedure below reported could be
(c(1:120),0.1,0.2,0.5,0.01,0.02,0.05,1)
thank you all of you for your attention
Stefano
function(y, trend1, trend2, seas, alphapar, beta, obsvar, maxiteration)
{
alpha <- alphapar
n <- length(y)
z <- vector("numeric", length = n)
for(timestep in 1:n) {
if(y[timestep] == 0) {
z[timestep] <- 0
}
else {
z[timestep] <- y[timestep]
}
}
# MATRICES THAT MUST BE MANUALLY INSERTED
# REGRESSION VECTOR F THAT BELONGS TO THE OBSERVATION EQUATION
FF <- matrix(data = c(1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), nrow =
14,
ncol = 1)
# THE CYCLICAL MATRIX G
G <- matrix(data = c(1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0,
0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0,
0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1,
0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0,
0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0,
0, 0, 0,
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0,
0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1,
0, 0, 0, 0,
0, 0, 0,
0, 0, 0, 0), nrow = 14, ncol = 14)
# THE VARIANCE-COVARIANCE MATRIX OF THE ERROR OF THE SYSTEM VECTOR:
# IT IS CONSTANT AND THE ROWS AND COLUMNS OF THE SEASONAL PART SUM TO ZERO
Vomega <- matrix(data = c(trend1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0,
trend2, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, seas, - seas/11, - seas/11,
- seas/11,
- seas/
11, - seas/11, - seas/11, - seas/11, - seas/11, - seas/11,
- seas/11,
- seas/
11, 0, 0, - seas/11, seas, - seas/11, - seas/11, - seas/11,
- seas/11,
- seas/
11, - seas/11, - seas/11, - seas/11, - seas/11, - seas/11,
0, 0, -
seas/11, -
seas/11, seas, - seas/11, - seas/11, - seas/11, - seas/11,
- seas/11,
- seas/
11, - seas/11, - seas/11, - seas/11, 0, 0, - seas/11, -
seas/11, -
seas/11,
seas, - seas/11, - seas/11, - seas/11, - seas/11, -
seas/11, -
seas/11, -
seas/11, - seas/11, 0, 0, - seas/11, - seas/11, - seas/11,
- seas/11,
seas, -
seas/11, - seas/11, - seas/11, - seas/11, - seas/11, -
seas/11, -
seas/11, 0,
0, - seas/11, - seas/11, - seas/11, - seas/11, - seas/11,
seas, -
seas/11, -
seas/11, - seas/11, - seas/11, - seas/11, - seas/11, 0, 0,
- seas/11,
- seas/
11, - seas/11, - seas/11, - seas/11, - seas/11, seas, -
seas/11, -
seas/11, -
seas/11, - seas/11, - seas/11, 0, 0, - seas/11, - seas/11,
- seas/11,
- seas/
11, - seas/11, - seas/11, - seas/11, seas, - seas/11, -
seas/11, -
seas/11, -
seas/11, 0, 0, - seas/11, - seas/11, - seas/11, - seas/11,
- seas/11,
- seas/
11, - seas/11, - seas/11, seas, - seas/11, - seas/11, -
seas/11, 0, 0,
- seas/
11, - seas/11, - seas/11, - seas/11, - seas/11, - seas/11,
- seas/11,
- seas/
11, - seas/11, seas, - seas/11, - seas/11, 0, 0, - seas/11,
- seas/11,
- seas/
11, - seas/11, - seas/11, - seas/11, - seas/11, - seas/11,
- seas/11,
- seas/
11, seas, - seas/11, 0, 0, - seas/11, - seas/11, - seas/11,
- seas/11,
- seas/
11, - seas/11, - seas/11, - seas/11, - seas/11, - seas/11,
- seas/11,
seas),
nrow = 14, ncol = 14)
inversecompoundsymmetric <- function(originalmatrix)
{
# par1 IS a
# par2 IS c
# par3 IS d
par1 <- originalmatrix[1, 1] - originalmatrix[1, 2]
par2 <- originalmatrix[1, 2]/par1
par3 <- par2/(1 + dim(originalmatrix)[[1]] * par2)
J <- matrix(1, nrow = dim(originalmatrix)[[1]], ncol =
dim(originalmatrix)[[1]])
inverseoriginalmatrix <- (diag(1, dim(originalmatrix)[[1]]) -
par3 *
J)/par1
result <- inverseoriginalmatrix
result
}
# matrixcompose CREATES A TWO-BLOCKS MATRIX
matrixcompose <- function(mat1, mat2)
{
dimtot <- dim(mat1)[[1]] + dim(mat2)[[1]]
dimmat1 <- dim(mat1)[[1]]
M <- matrix(0, nrow = dimtot, ncol = dimtot)
M[1:dimmat1, 1:dimmat1] <- mat1
M[(dimmat1 + 1):dimtot, (dimmat1 + 1):dimtot] <- mat2
result <- M
result
}
# THE UNKNOWN VARIABLE theta_t=1..n (OF DIMENSION (14*1))
thetasupport <- array(NA, dim = c(n, 14, 1))
theta <- array(NA, dim = c(n, 14, 1))
# E0 IS E[theta_1]
E0 <- array(NA, dim = c(n, 14, 1)) # V0 IS V[theta_1]
V0 <- array(NA, dim = c(n, 14, 14))
# E1 IS E[thetatilde_2]
E1 <- array(NA, dim = c(13, 1))
# V1 IS V[thetatilde_2]
V1 <- array(NA, dim = c(13, 13))
# CoV11 IS CoV[theta_t=1,thetatilde_t=2]
CoV11 <- array(NA, dim = c(14, 13))
# E4iniz IS E[theta_1|thetatilde_2]
E4iniz <- array(NA, dim = c(14, 1))
# V4iniz IS V[theta_1|thetatilde_2]
V4iniz <- array(NA, dim = c(14, 14))
# E5iniz IS E[psi_1|thetatilde_2]
E5iniz <- array(NA, dim = c(1, 1))
# V5 IS V[psi_t|thetatilde_2]
V5iniz <- array(NA, dim = c(1, 1))
# E6 IS E[thetastar_1|thetatilde_2]
E6iniz <- matrix(NA, nrow = 13, ncol = 1)
# V6 IS V[thetastar_1|thetatilde_2]
V6iniz <- array(NA, dim = c(13, 13))
# SETTING OF THE INITIAL CONDITIONS FOR ALL thetas, FROM TIME STEP 1 TO TIME
STEP n
for(timestep in 1:n) {
theta[timestep, 1, ] <- log((y[timestep] + 2), base = exp(1))
theta[timestep, (2:14), ] <- 0
}
# ITERATIONS
for(iterationumber in 1:maxiteration) {
# FIRST TIME-STEP
E0[1, 1:2, ] <- c(0.1, 0.01)
for(i in 3:14) {
E0[1, i, ] <- mean(y[1:12]) - y[i]
}
V0[1, , ] <- beta * Vomega
E1[, ] <- (G %*% E0[1, , ])[1:13]
V1[, ] <- (V0[1, , ] - V0[1, , ] %*% t(G) + Vomega)[1:13,
1:13]
CoV11 <- (V0[1, , ] %*% t(G))[1:14, 1:13]
E4iniz <- E0[1, , ] + CoV11 %*%
matrixcompose(inversecompoundsymmetric(V1[1:2, 1:
2]), inversecompoundsymmetric(V1[3:13, 3:13])) %*%
(theta[1:13] - E1)
V4iniz <- V0[1, , ] - CoV11 %*%
matrixcompose(inversecompoundsymmetric(V1[1:2, 1:
2]), inversecompoundsymmetric(V1[3:13, 3:13])) %*%
t(CoV11)
E5iniz <- E4iniz[1, 1] + E4iniz[3, 1]
V5iniz <- V4iniz[1, 1] + V4iniz[3, 3] + 2 * V4iniz[1, 3]
E6iniz <- E4iniz[(2:14), ]
V6iniz <- V4iniz[(2:14), (2:14)]
}
result <- list(E4iniz, E6iniz)
result
}
|