From 5b548cd6b56815c6bd922dc7bf35cbc88d042b54 Mon Sep 17 00:00:00 2001 From: Sebastian Krantz Date: Sun, 18 May 2025 02:29:09 +0200 Subject: [PATCH 1/2] Don't test on CRAN: runs too long. --- R/DFM.R | 4 ++-- man/DFM.Rd | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/R/DFM.R b/R/DFM.R index 0cf3f60..4e91ae3 100644 --- a/R/DFM.R +++ b/R/DFM.R @@ -138,7 +138,7 @@ #' BM14[, BM14_Models$freq == "M"] %<>% diff() #' BM14[, BM14_Models$freq == "Q"] %<>% diff(3) #' -#' +#' \donttest{ #' ### Small Model --------------------------------------- #' #' # IC for number of factors @@ -188,7 +188,7 @@ #' #' #' ### Large Model --------------------------------- -#' \donttest{ +#' #' # IC for number of factors #' IC_large <- ICr(BM14) #' plot(IC_large) diff --git a/man/DFM.Rd b/man/DFM.Rd index b27d344..c8bd046 100644 --- a/man/DFM.Rd +++ b/man/DFM.Rd @@ -147,7 +147,7 @@ BM14[, BM14_Models$log_trans] \%<>\% log() BM14[, BM14_Models$freq == "M"] \%<>\% diff() BM14[, BM14_Models$freq == "Q"] \%<>\% diff(3) - +\donttest{ ### Small Model --------------------------------------- # IC for number of factors @@ -197,7 +197,7 @@ plot(predict(dfm_medium), xlim = c(300, 370)) ### Large Model --------------------------------- -\donttest{ + # IC for number of factors IC_large <- ICr(BM14) plot(IC_large) From 0adcfedd6158e800df2a90c8eb2100a07ba023dd Mon Sep 17 00:00:00 2001 From: Sebastian Krantz Date: Tue, 17 Jun 2025 02:19:32 +0200 Subject: [PATCH 2/2] Add baseline idio MQ implementation (needs to be finished and verified). --- R/EMBM_MQ_idio.R | 133 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 133 insertions(+) create mode 100644 R/EMBM_MQ_idio.R diff --git a/R/EMBM_MQ_idio.R b/R/EMBM_MQ_idio.R new file mode 100644 index 0000000..1641c51 --- /dev/null +++ b/R/EMBM_MQ_idio.R @@ -0,0 +1,133 @@ +# Define the inputs of the function +# X: a matrix with n rows and TT columns representing the observed data at each time t. +# A: a transition matrix with size r*p+r*r. +# C: an observation matrix with size n*r. +# Q: a covariance matrix for the state equation disturbances with size r*p+r*r. +# R: a covariance matrix for the observation disturbances with size n*n. +# Z_0: the initial value of the state variable. +# V_0: the initial value of the variance-covariance matrix of the state variable. +# r: order of the state autoregression (AR) process. +# p: number of lags in the state AR process. +# i_idio: a vector of indicators specifying which variables are idiosyncratic noise variables.' + +# TODO: FINISH and VERIFY!! +# Z_0 = F_0; V_0 = P_0 +EMstepBMMQidio = function(X, A, C, Q, R, Z_0, V_0, XW0, NW, dgind, dnkron, dnkron_ind, r, p, R_mat, n, nq, sr, TT, rQi, rRi) { + + # Define dimensions of the input arguments + rC = 5L # ncol(Rcon) + pC = max(p, rC) + rp = r * p + rpC = r * pC + nm = n - nq + rC = rC * r # Note here replacing rC!! + + snm = seq_len(nm) + srp = seq_len(rp) + srpC = seq_len(rpC) + rpC1nq = (rpC+1L):(rpC+nq) + + # Run Kalman filter with current parameter estimates + kfs_res = SKFS(X, A, C, Q, R, Z_0, V_0, TRUE) + + # Extract values from the Kalman Filter output + Zsmooth = kfs_res$F_smooth + Vsmooth = kfs_res$P_smooth + VVsmooth = kfs_res$PPm_smooth + Zsmooth0 = kfs_res$F_smooth_0 + Vsmooth0 = kfs_res$P_smooth_0 + loglik = kfs_res$loglik + + # Normal Expectations for factors + tmp = rbind(Zsmooth0[srpC], Zsmooth[-TT, srpC, drop = FALSE]) + tmp2 = sum3(Vsmooth[srpC, srpC, -TT, drop = FALSE]) + EZZ = crossprod(Zsmooth[, srpC, drop = FALSE]) %+=% (tmp2 + Vsmooth[srpC, srpC, TT]) # E(Z'Z) + EZZ_BB = crossprod(tmp) %+=% (tmp2 + Vsmooth0[srpC, srpC]) # E(Z(-1)'Z(-1)) + EZZ_FB = crossprod(Zsmooth[, srpC, drop = FALSE], tmp) %+=% sum3(VVsmooth[srpC, srpC,, drop = FALSE]) # E(Z'Z(-1)) + + # Expectations for idiosyncratic errors + tmp = rbind(Zsmooth0[rpC1nq], Zsmooth[-TT, rpC1nq, drop = FALSE]) + tmp2 = sum3(Vsmooth[rpC1nq, rpC1nq, -TT, drop = FALSE]) + EZZ_u = diag(crossprod(Zsmooth[, rpC1nq, drop = FALSE])) + diag(tmp2 + Vsmooth[rpC1nq, rpC1nq, TT]) # E(Z'Z) + EZZ_BB_u = diag(crossprod(tmp)) + diag(tmp2 + Vsmooth0[rpC1nq, rpC1nq]) # E(Z(-1)'Z(-1)) + EZZ_FB_u = diag(diag(crossprod(Zsmooth[, rpC1nq, drop = FALSE], tmp)) + diag(rowSums(VVsmooth[rpC1nq, rpC1nq,, drop = FALSE], dims = 2L))) # E(Z'Z(-1)) + + # Update matrices A and Q + A_new = A + Q_new = Q + + # System matrices for factors + A_new[sr, srp] = EZZ_FB[sr, , drop = FALSE] %*% ainv(EZZ_BB) + if(rQi) { + Qsr = (EZZ[sr, sr] - tcrossprod(A_new[sr, srp, drop = FALSE], EZZ_FB[sr,, drop = FALSE])) / TT + Q_new[sr, sr] = if(rQi == 2L) Qsr else diag(diag(Qsr)) + } else Q_new[sr, sr] = diag(r) + + # System matrices for errors + A_new[rpC1nq, rpC1nq] = EZZ_FB_u / EZZ_BB_u + if(rRi) { + if(rRi == 2L) stop("Cannot estimate unrestricted observation covariance matrix together with AR(1) serial correlation") + Q_new[rpC1nq, rpC1nq] = (diag(EZZ_u) - A_new[rpC1nq, rpC1nq] * EZZ_FB_u) / TT + } else Q_new[rpC1nq, rpC1nq] = diag(nm) + + # Estimate matrix C using maximum likelihood approach + denom = numeric(nm*r^2) + nom = matrix(0, nm, r) + + for (t in seq_len(TT)) { + nmiss = as.double(NW[t, 1:nm]) + tmp = t(Zsmooth[t, sr]) + tmp2 = crossprod(tmp) + Vsmooth[sr, sr, t] + dim(tmp2) = NULL + denom %+=% tcrossprod(tmp2, nmiss) + nom %+=% (XW0[t, 1:nm] %*% tmp) + nom %-=% ((Zsmooth[t, rpC1nq] %*% tmp + Vsmooth[rpC1nq, sr, t]) * nmiss) + } + + dim(denom) = c(r, r, nm) + dnkron[dnkron_ind] = aperm.default(denom, c(1L, 3L, 2L)) + vec_C = solve.default(dnkron, unattrib(nom)) + C_new = C + C_new[1:nm, sr] = vec_C + + # Now updating the quarterly observation matrix C + for (i in (nm+1):n) { + denom = numeric(rC^2) + nom = numeric(rC) + i_idio_jQ = rpC + nm + 5*((i-nm)-1) + 1:5 + + for (t in 1:TT) { + if(NW[t, i]) { + denom %+=% (crossprod(Zsmooth[t, 1:rC, drop = FALSE]) + Vsmooth[1:rC, 1:rC, t]) + nom %+=% (XW0[t, i] * Zsmooth[t, 1:rC]) + nom %-=% (c(1,2,3,2,1) %*% (Zsmooth[t, i_idio_jQ] %*% t(Zsmooth[t, 1:rC]) + Vsmooth[i_idio_jQ, 1:rC, t])) + } + } + + dim(denom) = c(rC, rC) + denom_inv = ainv(denom) + C_i = denom_inv %*% nom + tmp = tcrossprod(denom_inv, R_mat) + C_i_constr = C_i - tmp %*% ainv(R_mat %*% tmp) %*% R_mat %*% C_i + C_new[i, 1:rC] = C_i_constr + + # Update AR parameters for quarterly idiosyncratic errors + V_0_new = V_0 + 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_new[rpC1nq[i-nm], rpC1nq[i-nm]] + Q_new[i_idio_jQ[1], i_idio_jQ[1]] = Q_new[rpC1nq[i-nm], rpC1nq[i-nm]] + } + + # Set initial conditions + V_0_new = V_0 + V_0_new[srpC, srpC] = Vsmooth0[srpC, srpC] + V_0_new[rpC1nq, rpC1nq] = diag(diag(Vsmooth0[rpC1nq, rpC1nq, drop = FALSE])) + + return(list(A = A_new, + C = C_new, + Q = Q_new, + R = R, # R stays fixed + F_0 = drop(Zsmooth0), + P_0 = V_0_new, + loglik = loglik)) +}