Nothing
#' @importFrom magic adiag
remNaNs_spline <-function(X,options){
TT <- dim(X)[1]
N <- dim(X)[2]
k <- options$k
indNaN <- is.na(X)
if(options$method == 1){ # replace all the missing values (método Giannone et al. 2008)
for (i in 1:N){
x = X[,i]
x[indNaN[,i]] = median(x,na.rm = T);
x_MA<-filter(x = c(rep(x[1],k),x,rep(x[length(x)],k)),filter = rep(1,2*k+1)/(2*k+1),sides = 1)
x_MA=x_MA[(2*k+1):length(x_MA)]
x[indNaN[,i]]=x_MA[indNaN[,i]]
X[,i]=x;
}
}else if(options$method == 2){ # replace missing values after removing leading and closing zeros
rem1 <- (rowSums(indNaN)>N*0.8)
nanLead <- which(rem1)
# nanEnd <- which(rem1[length(rem1):1])
# nanLE <- c(nanEnd,nanLead)
nanLE<-nanLead
X<-X[-nanLE,]
indNaN=is.na(X)
for (i in 1:N){
x = X[,i]
isnanx = is.na(x)
t1 = min(which(!isnanx))
t2 = max(which(!isnanx))
x1<-stats::spline(x[t1:t2],xout = 1:(t2-t1+1))
xx<-x1$y
x[t1:t2]<-x1$y
isnanx<-is.na(x)
x[isnanx] <- median(x,na.rm = T)
x_MA<-filter(x = c(rep(x[1],k),x,rep(x[length(x)],k)),filter = rep(1,2*k+1)/(2*k+1),sides = 1)
x_MA=x_MA[(2*k+1):length(x_MA)]
x[indNaN[,i]]=x_MA[indNaN[,i]]
X[,i]=x;
}
}else if(options$method == 3){ # only remove rows with leading and closing zeros
rem1 <- (rowSums(indNaN)==N)
nanLead <- which(rem1)
# nanEnd <- which(rem1[length(rem1):1])
# nanLE <- c(nanEnd,nanLead)
nanLE<-nanLead
if(length(nanLE) != 0){
X <- X[-nanLE,]
}
indNaN=is.na(X)
}else if(options$method == 4){ # remove rows with leading and closing zeros & replace missing values
rem1 <- (rowSums(indNaN)==N)
nanLead <- which(rem1)
# nanEnd <- which(rem1[length(rem1):1])
# nanLE <- c(nanEnd,nanLead)
nanLE<-nanLead
X<-X[-nanLE,]
indNaN=is.na(X)
for (i in 1:N){
x = X[,i]
isnanx = is.na(x)
t1 = min(which(!isnanx))
t2 = max(which(!isnanx))
x1<-stats::spline(x[t1:t2],xout = 1:(t2-t1+1))
xx<-x1$y
x[t1:t2]<-x1$y
isnanx<-is.na(x)
x[isnanx] <- median(x,na.rm = T)
x_MA<-filter(x = c(rep(x[1],k),x,rep(x[length(x)],k)),filter = rep(1,2*k+1)/(2*k+1),sides = 1)
x_MA=x_MA[(2*k+1):length(x_MA)]
x[indNaN[,i]]=x_MA[indNaN[,i]]
X[,i]=x;
}
}else if(options$method == 5){
indNaN=is.na(X)
for (i in 1:N){
x = X[,i]
isnanx = is.na(x)
t1 = min(which(!isnanx))
t2 = max(which(!isnanx))
x1<-stats::spline(x[t1:t2],xout = 1:(t2-t1+1))
xx<-x1$y
x[t1:t2]<-x1$y
isnanx<-is.na(x)
x[isnanx] <- median(x,na.rm = T)
x_MA<-filter(x = c(rep(x[1],k),x,rep(x[length(x)],k)),filter = rep(1,2*k+1)/(2*k+1),sides = 1)
x_MA=x_MA[(2*k+1):length(x_MA)]
x[indNaN[,i]]=x_MA[indNaN[,i]]
X[,i]=x;
}
}
return(list(X = X,indNaN=indNaN))
}
InitCond<-function(xNaN,r,p,blocks,optNaN,R_mat,q,nQ,i_idio){
x<-xNaN
Rcon<-R_mat
# library(magic)
pC = size(Rcon,2)
ppC = max(p,pC)
n_b = size(blocks,2)
res_remNaNs_spline <- remNaNs_spline(x,optNaN)
xBal<- res_remNaNs_spline$X
indNaN <- res_remNaNs_spline$indNaN
TT <- dim(xBal)[1] # time T of series
N <- dim(xBal)[2] # number of series
NM <- N-nQ # number of monthy frequency series
xNaN = xBal
for(i in 1:N){
xNaN[indNaN[,i],i] <- NA
}
# Initialize model coefficient output
C = {}
A = {}
Q = {}
initV = {}
res = xBal
resNaN = xNaN
# Set the first observations as NaNs: For quarterly-monthly aggreg. scheme
indNaN[1:pC-1,] <- T
for(i in 1:n_b){ # loop for each block
r_i<-r[i]
########################
# Observation equation #
########################
C_i = zeros(N,r_i*ppC) # Initialize state variable matrix helper
idx_i = find(blocks[,i]) # returns the series loaded in block i
# Note that the variables have been reshuffled so as to have all quarterly at the last columns of X
idx_iM = idx_i[idx_i<NM+1]; # index for monthly variables
idx_iQ = idx_i[idx_i>NM]; # index for quarterly variables
eig<-eigen(cov(res[,idx_iM]))
v<-eig$vectors[,1:r_i]
d<-eig$values[1:r_i]
C_i[idx_iM,1:r_i] = v
f = as.matrix(res[,idx_iM])%*%as.matrix(v)
for(kk in 0:(max(p+1,pC)-1)){
if(kk == 0){
FF<-f[(pC-kk):(dim(f)[1]-kk),]
}else{
FF <- cbind(FF,f[(pC-kk):(dim(f)[1]-kk),])
}
}
Rcon_i = kronecker(Rcon,eye(r_i))
q_i = kronecker(q,zeros(r_i,1));
ff = FF[,1:(r_i*pC)]
for(j in idx_iQ){ # Coeficiente "loadings" de Variáveis trimestrais
xx_j = resNaN[pC:dim(resNaN)[1],j]
if(sum(!is.na(xx_j)) < size(ff,2)+2){
xx_j = res[pC:dim(res)[1],j]
}
ff_j = ff[!is.na(xx_j),]
xx_j = xx_j[!is.na(xx_j)]
iff_j = solve(t(ff_j)%*%ff_j)
Cc = iff_j%*%t(ff_j)%*%xx_j
Cc = Cc - iff_j%*%t(Rcon_i)%*%solve(Rcon_i%*%iff_j%*%t(Rcon_i))%*%(Rcon_i%*%Cc-q_i);
C_i[j,1:(pC*r_i)] <- t(Cc)
}
ff = rbind(zeros(pC-1,pC*r_i),ff)
res = res - ff%*%t(C_i)
resNaN = res
for(i_aux in 1:dim(indNaN)[2]){
resNaN[indNaN[,i_aux],i_aux] <- NA
}
C <- cbind(C,C_i)
#######################
# Transition Equation #
#######################
z <- FF[,1:r_i]
Z <- FF[,(r_i+1):(r_i*(p+1))]
A_i = t(zeros(r_i*ppC,r_i*ppC))
A_temp = solve(t(Z)%*%Z)%*%t(Z)%*%z
A_i[1:r_i,1:(r_i*p)] <- t(A_temp)
A_i[(r_i+1):dim(A_i)[1],1:(r_i*(ppC-1))] <- eye(r_i*(ppC-1))
##########################
Q_i = zeros(ppC*r_i,ppC*r_i)
e = z - Z%*%A_temp # VAR residuals
Q_i[1:r_i,1:r_i] = cov(e); # VAR covariance matrix
initV_i = matrix(solve(eye((r_i*ppC)^2)-kronecker(A_i,A_i))%*%c(Q_i),r_i*ppC,r_i*ppC);
if(is.null(A)){
A<-A_i
}else{
A <- magic::adiag(A,A_i)
}
if(is.null(Q)){
Q<-Q_i
}else{
Q <- magic::adiag(Q,Q_i)
}
if(is.null(initV)){
initV<-initV_i
}else{
initV <- magic::adiag(initV,initV_i)
}
}
R = diag(apply(resNaN, 2, stats::var, na.rm = T))
eyeN = eye(N)
eyeN<-eyeN[,i_idio]
# Initial conditions
C <- cbind(C,eyeN)
ii_idio = find(i_idio)
n_idio = length(ii_idio)
B = zeros(n_idio)
S = zeros(n_idio)
BM <- zeros(n_idio)
SM <- zeros(n_idio)
# Loop for monthly variables
for (i in 1:n_idio){
# Set observation equation residual covariance matrix diagonal
R[ii_idio[i],ii_idio[i]] <- 1e-04
# Subsetting series residuals for series i
res_i <- resNaN[,ii_idio[i]]
# number of leading zeros
# leadZero = max( find( t(1:TT) == cumsum(is.na(res_i)) ) )
# endZero = max( find( TT:1 == cumsum(is.na(res_i[length(res_i):1])) ) );
# res_i<-res_i[(leadZero+1):(length(res_i)-endZero)]
res_i<-res_i[!is.na(res_i)]
# Linear regression: AR 1 process for monthly series residuals
BM[i,i] = solve(t(res_i[1:(length(res_i)-1)])%*%res_i[1:(length(res_i)-1)])%*%t(res_i[1:(length(res_i)-1)])%*%res_i[2:length(res_i)]
SM[i,i] = stats::var(res_i[2:length(res_i)]-res_i[1:(length(res_i)-1)]*BM[i,i])
}
# blocks for covariance matrices
initViM = diag(1/diag(eye(size(BM,1))-BM^2))%*%SM;
if(!nQ==0){
C<-cbind(C,rbind(zeros(NM,5*nQ),t(kronecker(eye(nQ),c(1,2,3,2,1)))))
Rdiag<-diag(R)
sig_e <- Rdiag[(NM+1):N]/19
Rdiag[(NM+1):N] <- 1e-04
R <- diag(Rdiag) # Covariance for obs matrix residuals
# for BQ, SQ
rho0 <- 0.1
temp <- zeros(5)
temp[1,1] <- 1
# BQ and SQ
BQ <- kronecker(eye(nQ),rbind(cbind(rho0,zeros(1,4)),cbind(eye(4),zeros(4,1))))
if(length(sig_e)>1){
SQ = kronecker(diag((1-rho0^2)*sig_e),temp)
}else{
SQ = kronecker((1-rho0^2)*sig_e,temp)
}
initViQ = matrix(solve(eye((5*nQ)^2)-kronecker(BQ,BQ))%*%c(SQ),5*nQ,5*nQ)
}else{
BQ <- NULL
SQ <- NULL
initViQ <- NULL
}
A1<-magic::adiag(A,BM,BQ)
Q1<-magic::adiag(Q, SM, SQ)
A<-A1 # Observation matrix
Q<-Q1 # Residual covariance matrix
# Initial conditions
initZ = zeros(size(A,1),1); # States
initV = magic::adiag(initV, initViM, initViQ) # Covariance of states
return(list(A = A, C = C, Q = Q, R = R, initZ = initZ, initV = initV))
}
EM_DFM_SS_block_idioQARMA_restrMQ<-function(X,Par){
# library(matlab)
thresh = 1e-4
r = Par$r
p = Par$p
max_iter = Par$max_iter
i_idio = Par$i_idio
R_mat = Par$Rconstr
q = Par$q
nQ = Par$nQ
blocks = Par$blocks
### Prepara??o dos dados
TT <- dim(X)[1]
N <- dim(X)[2]
### Standardise X
Mx = colMeans(X,na.rm=T)
Wx = sapply(1:N,function(x) sd(X[,x],na.rm = T))
xNaN <- (X-kronecker(t(Mx),rep(1,TT)))/kronecker(t(Wx),rep(1,TT))
### Initial conditions
# Removing missing values (for initial estimator)
optNaN<-list()
optNaN$method = 2; # Remove leading and closing zeros
optNaN$k = 3;
res_InitCond<-InitCond(xNaN,r,p,blocks,optNaN,R_mat,q,nQ,i_idio);
A<-res_InitCond$A
C<-res_InitCond$C
Q<-res_InitCond$Q
R<-res_InitCond$R
Z_0<-res_InitCond$initZ
V_0<-res_InitCond$initV
# some auxiliary variables for the iterations
previous_loglik = -Inf
num_iter = 0
LL = -Inf
converged = F
# y for the estimation is WITH missing data
y = t(xNaN)
#--------------------------------------------------------------------------
#THE EM LOOP
#--------------------------------------------------------------------------
#The model can be written as
#y = C*Z + e;
#Z = A*Z(-1) + v
#where y is NxT, Z is (pr)xT, etc
#remove the leading and ending nans for the estimation
optNaN$method = 3
y_est <- remNaNs_spline(xNaN,optNaN)
y_est_indNaN<-t(y_est$indNaN)
y_est<-t(y_est$X)
num_iter = 0
loglik_aux = -Inf
while ((num_iter < max_iter) & !converged){
res_EMstep = EMstep(y_est, A, C, Q, R, Z_0, V_0, r,p,R_mat,q,nQ,i_idio,blocks)
# res_EMstep <- list(C_new, R_new, A_new, Q_new, Z_0, V_0, loglik)
C = res_EMstep$C_new;
R = res_EMstep$R_new;
A = res_EMstep$A_new;
Q = res_EMstep$Q_new;
Z_0<-res_EMstep$Z_0
V_0<-res_EMstep$V_0
loglik<-res_EMstep$loglik
# updating the user on what is going on
if((num_iter%% 5)==0&&num_iter>0){message(num_iter,"th iteration: \nThe loglikelihood went from ",round(loglik_aux,4)," to ",round(loglik,4))}
loglik_aux <- loglik
# Checking convergence
if (num_iter>2){
res_em_converged = em_converged(loglik, previous_loglik, thresh,1)
converged<-res_em_converged$converged
decreasse<-res_em_converged$decrease
}
LL <- cbind(LL, loglik)
previous_loglik <- loglik
num_iter <- num_iter + 1
}
# final run of the Kalman filter
res_runKF = runKF(y, A, C, Q, R, Z_0, V_0)
Zsmooth<-t(res_runKF$xsmooth)
x_sm <- Zsmooth[2:dim(Zsmooth)[1],]%*%t(C)
Res<-list()
Res$x_sm <- x_sm
Res$X_sm <- kronecker(t(Wx),rep(1,TT))*x_sm + kronecker(t(Mx),rep(1,TT))
Res$FF <- Zsmooth[2:dim(Zsmooth)[1],]
#--------------------------------------------------------------------------
# Loading the structure with the results
#--------------------------------------------------------------------------
Res$C <- C;
Res$R <- R;
Res$A <- A;
Res$Q <- Q;
Res$Mx <- Mx;
Res$Wx <- Wx;
Res$Z_0 <- Z_0;
Res$V_0 <- V_0;
Res$r <- r;
Res$p <- p;
return(Res)
}
EMstep <- function(y = NULL, A = NULL, C = NULL, Q = NULL, R = NULL, Z_0 = NULL, V_0 = NULL,
r = NULL, p = NULL, R_mat = NULL, q = NULL, nQ = NULL, i_idio = NULL, blocks = NULL){
n <- size(y,1)
TT <- size(y,2)
nM <- n - nQ
pC <- size(R_mat,2)
ppC <- max(p,pC)
n_b <- size(blocks,2)
# Compute the (expected) sufficient statistics for a single Kalman filter sequence.
#Running the Kalman filter with the current estimates of the parameters
res_runKF = runKF(y, A, C, Q, R, Z_0, V_0);
Zsmooth<-res_runKF$xsmooth
Vsmooth<-res_runKF$Vsmooth
VVsmooth<-res_runKF$VVsmooth
loglik<-res_runKF$loglik
A_new <- A
Q_new <- Q
V_0_new <- V_0
for(i in 1:n_b){
r_i <- r[i]
rp <- r_i*p
if(i == 1){
rp1 <- 0*ppC
}else{
rp1 <- sum(r[1:(i-1)])*ppC
}
A_i <- A[(rp1+1):(rp1+r_i*ppC), (rp1+1):(rp1+r_i*ppC)]
Q_i <- Q[(rp1+1):(rp1+r_i*ppC), (rp1+1):(rp1+r_i*ppC)]
if(rp==1){
EZZ <- t(Zsmooth[(rp1+1):(rp1+rp),2:ncol(Zsmooth)]) %*% Zsmooth[(rp1+1):(rp1+rp),2:ncol(Zsmooth)] +
sum(Vsmooth[(rp1+1):(rp1+rp),(rp1+1):(rp1+rp),2:dim(Vsmooth)[3]]) # E(Z'Z)
EZZ_BB <- t(Zsmooth[(rp1+1):(rp1+rp),1:(ncol(Zsmooth)-1)]) %*% Zsmooth[(rp1+1):(rp1+rp),1:(ncol(Zsmooth)-1)] +
sum(Vsmooth[(rp1+1):(rp1+rp),(rp1+1):(rp1+rp),1:(dim(Vsmooth)[3]-1)]) #E(Z(-1)'Z_(-1))
EZZ_FB <- t(Zsmooth[(rp1+1):(rp1+rp),2:ncol(Zsmooth)]) %*% Zsmooth[(rp1+1):(rp1+rp),1:(ncol(Zsmooth)-1)] +
sum(VVsmooth[(rp1+1):(rp1+rp),(rp1+1):(rp1+rp),]) #E(Z'Z_(-1))
}else{
EZZ <- (Zsmooth[(rp1+1):(rp1+rp),2:ncol(Zsmooth)]) %*% t(Zsmooth[(rp1+1):(rp1+rp),2:ncol(Zsmooth)]) +
apply(Vsmooth[(rp1+1):(rp1+rp),(rp1+1):(rp1+rp),2:dim(Vsmooth)[3]],c(1,2),sum) # E(Z'Z)
EZZ_BB <- (Zsmooth[(rp1+1):(rp1+rp),1:(ncol(Zsmooth)-1)]) %*% t(Zsmooth[(rp1+1):(rp1+rp),1:(ncol(Zsmooth)-1)]) +
apply(Vsmooth[(rp1+1):(rp1+rp),(rp1+1):(rp1+rp),1:(dim(Vsmooth)[3]-1)],c(1,2),sum) #E(Z(-1)'Z_(-1))
EZZ_FB <- (Zsmooth[(rp1+1):(rp1+rp),2:ncol(Zsmooth)]) %*% t(Zsmooth[(rp1+1):(rp1+rp),1:(ncol(Zsmooth)-1)]) +
apply(VVsmooth[(rp1+1):(rp1+rp),(rp1+1):(rp1+rp),],c(1,2),sum) #E(Z'Z_(-1))
}
A_i[1:r_i,1:rp] <- EZZ_FB[1:r_i,1:rp] %*% solve(EZZ_BB[1:rp,1:rp])
Q_i[1:r_i,1:r_i] <- (EZZ[1:r_i,1:r_i] - A_i[1:r_i,1:rp] %*% t(matrix(EZZ_FB[1:r_i,1:rp],r_i,rp))) / TT
A_new[(rp1+1):(rp1+r_i*ppC),(rp1+1):(rp1+r_i*ppC)] <- A_i
Q_new[(rp1+1):(rp1+r_i*ppC),(rp1+1):(rp1+r_i*ppC)] <- Q_i;
V_0_new[(rp1+1):(rp1+r_i*ppC),(rp1+1):(rp1+r_i*ppC)] <- Vsmooth[(rp1+1):(rp1+r_i*ppC),(rp1+1):(rp1+r_i*ppC),1]
}
rp1 <- sum(sum(r))*ppC
niM <- sum(i_idio[1:nM])
# idiosyncratic
EZZ <- diag(diag(Zsmooth[(rp1+1):nrow(Zsmooth),2:ncol(Zsmooth)] %*% t(Zsmooth[(rp1+1):nrow(Zsmooth),2:ncol(Zsmooth)]))) + diag(diag(apply(Vsmooth[(rp1+1):dim(Vsmooth)[1],(rp1+1):dim(Vsmooth)[2],2:dim(Vsmooth)[3]], MARGIN = 1:2, FUN = sum))) #E(Z'Z)
EZZ_BB <- diag(diag(Zsmooth[(rp1+1):nrow(Zsmooth),1:(ncol(Zsmooth)-1)] %*% t(Zsmooth[(rp1+1):nrow(Zsmooth),1:(ncol(Zsmooth)-1)]))) + diag(diag(apply(Vsmooth[(rp1+1):dim(Vsmooth)[1],(rp1+1):dim(Vsmooth)[2],1:(dim(Vsmooth)[3]-1)], MARGIN = 1:2, FUN = sum))) #E(Z(-1)'Z_(-1))
EZZ_FB <- diag(diag(Zsmooth[(rp1+1):nrow(Zsmooth),2:ncol(Zsmooth)] %*% t(Zsmooth[(rp1+1):nrow(Zsmooth),1:(ncol(Zsmooth)-1)]))) + diag(diag(apply(VVsmooth[(rp1+1):dim(VVsmooth)[1],(rp1+1):dim(VVsmooth)[2],], MARGIN = 1:2, FUN = sum))) #E(Z'Z_(-1))
A_i <- EZZ_FB %*% diag(1/diag(EZZ_BB))
Q_i <- (EZZ - A_i %*% t(EZZ_FB)) / TT
A_new[(rp1+1):(rp1+niM),(rp1+1):(rp1+niM)] <- A_i[1:niM,1:niM]
Q_new[(rp1+1):(rp1+niM),(rp1+1):(rp1+niM)] <- Q_i[1:niM,1:niM]
V_0_new[(rp1+1):(rp1+niM),(rp1+1):(rp1+niM)] <- diag(diag(Vsmooth[(rp1+1):(rp1+niM),(rp1+1):(rp1+niM),1]))
Z_0 <- Zsmooth[,1] #zeros(size(Zsmooth,1),1); #
# nanY <- is.nan(y)
nanY<-is.na(y)
y[nanY] <- 0
# LOADINGS
C_new <- C
# Blocks
bl <- unique(blocks)
n_bl <- size(bl,1)
bl_idxM <- NULL
bl_idxQ <- NULL
R_con <- NULL
q_con <- NULL
for(i in 1:n_b){
bl_idxQ <- cbind(bl_idxQ, repmat(bl[,i],1,r[i]*ppC))
bl_idxM <- cbind(bl_idxM, repmat(bl[,i],1,r[i]), zeros(n_bl,r[i]*(ppC-1)))
if(i == 1){
R_con <- kronecker(R_mat,eye(r[i]))
}else{
R_con <- as.matrix(Matrix::bdiag(R_con, kronecker(R_mat,eye(r[i]))))
}
q_con <- rbind(q_con, zeros(r[i]*size(R_mat,1),1))
}
bl_idxM <- bl_idxM == 1
bl_idxQ <- bl_idxQ == 1
#idio
i_idio_M <- i_idio[1:nM]
n_idio_M <- length(find(i_idio_M))
c_i_idio <- cumsum(i_idio)
for(i in 1:n_bl){
bl_i <- bl[i,]
rs <- sum(r[bl_i == 1])
idx_i <- NULL
for(k in 1:nrow(blocks)){
idx_i[k] <- sum(blocks[k,] == bl_i) == size(blocks,2)
}
idx_i <- find(idx_i)
# MONTHLY
idx_iM <- idx_i[idx_i < (c(nM) + 1)]
n_i <- length(idx_iM)
if(n_i==0){
denom <- NULL
nom <- NULL
} else {
denom <- zeros(n_i*rs,n_i*rs)
nom <- zeros(n_i,rs)
i_idio_i <- i_idio_M[idx_iM] == 1
i_idio_ii <- c_i_idio[idx_iM]
i_idio_ii <- i_idio_ii[i_idio_i]
for(t in 1:TT){
nanYt <- diag(!nanY[idx_iM,t])
nn <- sum(bl_idxM[i,])
denom <- denom + kronecker(Zsmooth[bl_idxM[i,],t+1][1:nn] %*% t(Zsmooth[bl_idxM[i,],t+1][1:nn]) + Vsmooth[bl_idxM[i,],bl_idxM[i,],t+1][1:nn,1:nn], nanYt)
nom <- nom + y[idx_iM,t] %*% t(Zsmooth[bl_idxM[i,],t+1][1:nn]) - nanYt[,i_idio_i] %*% (Zsmooth[rp1+i_idio_ii,t+1] %*% t(Zsmooth[bl_idxM[i,],t+1][1:nn]) + Vsmooth[rp1+i_idio_ii,bl_idxM[i,],t+1][,1:nn])
}
vec_C <- solve(denom) %*% c(nom)
C_new[idx_iM,bl_idxM[i,]][,1:nn] <- matrix(vec_C,n_i,rs)
}
# QUARTERLY
idx_iQ <- idx_i[idx_i>c(nM)]
rps <- rs*ppC
R_con_i <- R_con[,bl_idxQ[i,]]
q_con_i <- q_con
no_c <- !(rowSums(R_con_i == 0) == ncol(R_con_i))
R_con_i <- R_con_i[no_c,]
q_con_i <- q_con_i[no_c,]
for(j in idx_iQ){
denom <- zeros(rps,rps)
nom <- zeros(1,rps)
idx_jQ <- j-c(nM)
if(i != 1){
i_idio_jQ <- (rp1+n_idio_M+5*(idx_jQ-1)+1):(rp1+n_idio_M+5*idx_jQ)
V_0_new[i_idio_jQ,i_idio_jQ] <- Vsmooth[i_idio_jQ,i_idio_jQ,1]
A_new[i_idio_jQ[1],i_idio_jQ[1]] <- A_i[i_idio_jQ[1]-rp1,i_idio_jQ[1]-rp1]
Q_new[i_idio_jQ[1],i_idio_jQ[1]] <- Q_i[i_idio_jQ[1]-rp1,i_idio_jQ[1]-rp1]
for(t in 1:TT){
nanYt <- as.vector(!nanY[j,t])*1
nn2 <- sum(bl_idxQ[i,])
denom <- denom + kronecker(Zsmooth[bl_idxQ[i,],t+1][1:nn2] %*% t(Zsmooth[bl_idxQ[i,],t+1][1:nn2]) + Vsmooth[bl_idxQ[i,],bl_idxQ[i,],t+1][1:nn2,1:nn2],nanYt)
nom <- nom + y[j,t] %*% t(Zsmooth[bl_idxQ[i,],t+1][1:nn2])
nom <- nom - nanYt %*% (matrix(c(1,2,3,2,1), nrow = 1) %*% Zsmooth[i_idio_jQ,t+1] %*% t(Zsmooth[bl_idxQ[i,],t+1][1:nn2]) +
matrix(c(1,2,3,2,1), nrow = 1) %*% Vsmooth[i_idio_jQ,bl_idxQ[i,],t+1][,1:nn2])
}
C_i <- solve(denom) %*% t(nom)
C_i_constr <- C_i - solve(denom) %*% t(R_con_i) %*% solve(R_con_i %*% solve(denom) %*% t(R_con_i)) %*% (R_con_i %*% C_i - q_con_i)
nn3 <- sum(bl_idxQ[i,])
C_new[j,bl_idxQ[i,]][1:nn3] <- C_i_constr
}
}
}
R_new <- zeros(n,n)
for(t in 1:TT){
nanYt <- diag(!nanY[,t])*1 == 1
R_new <- R_new + (y[,t] - nanYt %*% C_new %*% Zsmooth[,t+1]) %*% t(y[,t] - nanYt %*% C_new %*% Zsmooth[,t+1]) + nanYt %*% C_new %*% Vsmooth[,,t+1] %*% t(C_new) %*% nanYt + (eye(n)-nanYt) %*% R %*% (eye(n)-nanYt)
}
R_new <- R_new/TT
RR <- diag(R_new) #RR(RR<1e-2) = 1e-2;
RR[i_idio_M] <- 1e-04
if(nM<length(RR)){
RR[(nM+1):length(RR)] <- 1e-04
}
R_new <- diag(RR)
if(!is.matrix(Z_0)){
Z_0<-matrix(Z_0,length(Z_0),1)
}
# output
return(list(C_new = C_new, R_new = R_new, A_new = A_new, Q_new = Q_new, Z_0 = Z_0, V_0 = V_0, loglik = loglik))
}
FIS <- function(Y,Z,R,TT,Q,S){
# library(corpcor)
# %______________________________________________________________________
# % Fixed intervall smoother (see Harvey, 1989, p. 154)
# % FIS returns the smoothed state vector AmT and its covar matrix PmT
# % Use this in conjunction with function SKF
# %______________________________________________________________________
# % INPUT
# % Y Data (nobs x n)
# % S Estimates from Kalman filter SKF
# % S.Am : Estimates a_t|t-1 (nobs x m)
# % S.Pm : P_t|t-1 = Cov(a_t|t-1) (nobs x m x m)
# % S.AmU : Estimates a_t|t (nobs x m)
# % S.PmU : P_t|t = Cov(a_t|t) (nobs x m x m)
# % OUTPUT
# % S Smoothed estimates added to above
# % S.AmT : Estimates a_t|T (nobs x m)
# % S.PmT : P_t|T = Cov(a_t|T) (nobs x m x m)
# % S.PmT_1 : Cov(a_ta_t-1|T)
# % where m is the dim of state vector and t = 1 ...T is time
m<-dim(S$Am)[1]
nobs<-dim(S$Am)[2]
S$AmT = zeros(m,nobs+1)
S$PmT = array(0,c(m,m,nobs+1))
S$AmT[,nobs+1] <- S$AmU[,nobs+1]
S$PmT[,,nobs+1] <- S$PmU[,,nobs+1]
S$PmT_1<-array(0,c(m,m,nobs))
S$PmT_1[,,nobs] <- (eye(m)-S$KZ)%*%TT%*%S$PmU[,,nobs]
pinv<-corpcor::pseudoinverse(S$Pm[,,nobs])
J_2 <- S$PmU[,,nobs]%*%t(TT)%*%pinv
for (t in nobs:1){
PmU <- S$PmU[,,t]
Pm1 <- S$Pm[,,t]
P_T <- S$PmT[,,t+1]
P_T1 <- S$PmT_1[,,t]
J_1 <- J_2
S$AmT[,t] <- S$AmU[,t] + J_1%*%(S$AmT[,t+1] - TT%*%S$AmU[,t])
S$PmT[,,t] <- PmU + J_1%*%(P_T - Pm1)%*%t(J_1)
if(t>1){
pinv<-corpcor::pseudoinverse(S$Pm[,,t-1])
J_2 <- S$PmU[,,t-1]%*%t(TT)%*%pinv
S$PmT_1[,,t-1] = PmU%*%t(J_2)+J_1%*%(P_T1-TT%*%PmU)%*%t(J_2)
}
}
return(S)
}
SKF <-function(Y,Z,R,TT,Q,A_0,P_0){
#Y = y; Z = C; TT = A; A_0 = x_0; P_0 = Sig_0
# %______________________________________________________________________
# % Kalman filter for stationary systems with time-varying system matrices
# % and missing data.
# %
# % The model is y_t = Z * a_t + eps_t
# % a_t+1 = TT * a_t + u_t
# %
# %______________________________________________________________________
# % INPUT
# % Y Data (nobs x n)
# % OUTPUT
# % S.Am Predicted state vector A_t|t-1 (nobs x m)
# % S.AmU Filtered state vector A_t|t (nobs+1 x m)
# % S.Pm Predicted covariance of A_t|t-1 (nobs x m x m)
# % S.PmU Filtered covariance of A_t|t (nobs+1 x m x m)
# % S.loglik Value of likelihood function
#
# % Output structure & dimensions
n <- dim(Z)[1]
m <- dim(Z)[2]
nobs <- size(Y,2)
S<-list()
S$Am <- array(NA,c(m,nobs))
S$Pm <- array(NA,c(m,m,nobs))
S$AmU <- array(NA,c(m,nobs+1))
S$PmU <- array(NA,c(m,m,nobs+1))
S$loglik <- 0
# ______________________________________________________________________
Au <- A_0; # A_0|0;
Pu <- P_0; # P_0|0
S$AmU[,1] = Au;
S$PmU[,,1] = Pu;
for(t in 1:nobs){
# t
# A = A_t|t-1 & P = P_t|t-1
A <- TT%*%Au;
P <- TT%*%Pu%*%t(TT) + Q;
P <- 0.5*(P+t(P))
# handling the missing data
res_MissData <- MissData(Y[,t],Z,R)
y_t <- res_MissData$y
Z_t <- res_MissData$C
R_t <- res_MissData$R
L_t <- res_MissData$L
# if(is.null(y_t)){
if(sum(is.na(y_t))==length(y_t)){
Au <- A
Pu <- P
} else {
if(!is.matrix(Z_t)){
Z_t<-t(as.matrix(Z_t))
}
PZ <- P%*%t(Z_t)
iF <- solve(Z_t%*%PZ + R_t)
PZF <- PZ%*%iF
V <- y_t - Z_t%*%A
Au <- A + PZF%*%V
Pu <- P - PZF%*%t(PZ)
Pu <- 0.5*(Pu+t(Pu))
S$loglik <- S$loglik + 0.5*(log(det(iF)) - t(V)%*%iF%*%V)
}
S$Am[,t] <- A
S$Pm[,,t] <- P
# Au = A_t|t & Pu = P_t|t
S$AmU[,t+1] <- Au
S$PmU[,,t+1] <- Pu
} # t
if(sum(is.na(y_t))==length(y_t)){
S$KZ <- zeros(m,m)
}else{
S$KZ <- PZF%*%Z_t
}
return(S)
}
MissData <- function(y,C,R){
# library(matlab)
# ______________________________________________________________________
# PROC missdata
# PURPOSE: eliminates the rows in y & matrices Z, G that correspond to
# missing data (NaN) in y
# INPUT y vector of observations at time t (n x 1 )
# S KF system matrices (structure)
# must contain Z & G
# OUTPUT y vector of observations (reduced) (# x 1)
# Z G KF system matrices (reduced) (# x ?)
# L To restore standard dimensions (n x #)
# where # is the nr of available data in y
# ______________________________________________________________________
# [y,C,R,L]
ix <- !is.na(y)
e <- eye(length(y))
L <- e[,ix]
y <- y[ix]
C <- C[ix,]
R <- R[ix,ix]
return(list(y=y,C=C,R=R,L=L))
}
# %%% Replication files for:
# %%% ""Nowcasting", 2010, (by Marta Banbura, Domenico Giannone and Lucrezia Reichlin),
# %%% in Michael P. Clements and David F. Hendry, editors, Oxford Handbook on Economic Forecasting.
# %%%
# %%% The software can be freely used in applications.
# %%% Users are kindly requested to add acknowledgements to published work and
# %%% to cite the above reference in any resulting publications
# %--------------------------------------------------------------------------
# % KALMAN FILTER
# %--------------------------------------------------------------------------
runKF <- function(y, A, C, Q, R, x_0, Sig_0){
# x_0 = Z_0; Sig_0 = V_0
S <- SKF(y,C,R,A,Q, x_0, Sig_0);
S <- FIS(y,C,R,A,Q,S);
xsmooth <- S$AmT;
Vsmooth <- S$PmT;
VVsmooth <- S$PmT_1;
loglik <- S$loglik;
return(list(xsmooth = xsmooth,Vsmooth = Vsmooth,VVsmooth = VVsmooth,loglik = loglik))
}
em_converged <- function(loglik = NULL, previous_loglik = NULL, threshold = NULL, check_increased = NULL){
# EM_CONVERGED Has EM converged?
# [converged, decrease] = em_converged(loglik, previous_loglik, threshold)
#
# We have converged if the slope of the log-likelihood function falls below 'threshold',
# i.e., |f(t) - f(t-1)| / avg < threshold,
# where avg = (|f(t)| + |f(t-1)|)/2 and f(t) is log lik at iteration t.
# 'threshold' defaults to 1e-4.
#
# This stopping criterion is from Numerical Recipes in C p423
#
# If we are doing MAP estimation (using priors), the likelihood can decrase,
# even though the mode of the posterior is increasing.
nargin <- 4 - sum(c(is.null(loglik), is.null(previous_loglik), is.null(threshold), is.null(check_increased)))
if(nargin < 3){threshold <- 1e-4}
if(nargin < 4){check_increased <- 1}
converged <- 0
decrease <- 0
if(!is.null(check_increased)){
if(loglik - previous_loglik < -1e-3){ # allow for a little imprecision
message(paste("******likelihood decreased from",round(previous_loglik,4),"to",round(loglik,4)))
decrease <- 1
}
}
delta_loglik <- abs(loglik - previous_loglik)
avg_loglik <- (abs(loglik) + abs(previous_loglik) + 2.2204e-16)/2
if((delta_loglik / avg_loglik) < threshold){converged <- 1}
# output
list(converged = converged, decrease = decrease)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.