No update for the Initial State Covariance Matrix V_0

[victorbluhu]

I believe that in line 243 of the dfm.m file the EMstep function should return V_0_new and not V_0. Due to this typo, the initial state covariance matrix is not updated in the algorithm. It does not seem to affect estimations by a lot, but we are leaving loglikelihood on the table.

[MajesticKhan]

HI @victorbluhu ,

This is not a typo. According to the function's documentation, this is expected

Nowcasting/functions/dfm.m

Lines 244 to 292 in 19f365c

	%EMstep Applies EM algorithm for parameter reestimation
	%
	% Syntax:
	% [C_new, R_new, A_new, Q_new, Z_0, V_0, loglik]
	% = EMstep(y, A, C, Q, R, Z_0, V_0, r, p, R_mat, q, nQ, i_idio, blocks)
	%
	% Description:
	% EMstep reestimates parameters based on the Estimation Maximization (EM)
	% algorithm. This is a two-step procedure:
	% (1) E-step: the expectation of the log-likelihood is calculated using
	% previous parameter estimates.
	% (2) M-step: Parameters are re-estimated through the maximisation of
	% the log-likelihood (maximize result from (1)).
	%
	% See "Maximum likelihood estimation of factor models on data sets with
	% arbitrary pattern of missing data" for details about parameter
	% derivation (Banbura & Modugno, 2010). This procedure is in much the
	% same spirit.
	%
	% Input:
	% y: Series data
	% A: Transition matrix
	% C: Observation matrix
	% Q: Covariance for transition equation residuals
	% R: Covariance for observation matrix residuals
	% Z_0: Initial values of factors
	% V_0: Initial value of factor covariance matrix
	% r: Number of common factors for each block (e.g. vector [1 1 1 1])
	% p: Number of lags in transition equation
	% R_mat: Estimation structure for quarterly variables (i.e. "tent")
	% q: Constraints on loadings
	% nQ: Number of quarterly series
	% i_idio: Indices for monthly variables
	% blocks: Block structure for each series (i.e. for a series, the structure
	% [1 0 0 1] indicates loadings on the first and fourth factors)
	%
	% Output:
	% C_new: Updated observation matrix
	% R_new: Updated covariance matrix for residuals of observation matrix
	% A_new: Updated transition matrix
	% Q_new: Updated covariance matrix for residuals for transition matrix
	% Z_0: Initial value of state
	% V_0: Initial value of covariance matrix
	% loglik: Log likelihood
	%
	% References:
	% "Maximum likelihood estimation of factor models on data sets with
	% arbitrary pattern of missing data" by Banbura & Modugno (2010).
	% Abbreviated as BM2010

[victorbluhu]

Hi @MajesticKhan , sorry for the very late response.

Even though the code states "initial value of ____" (state / covariance matrix), they do update the Z_0 value with the smoothed post-Kalman Filter values Zsmooth(:,1) (as you may see in lines 312 (KF filter-smoother) and 403 (Z_0 update)). The Z_0 in the output is different from the Z_0 in the input.

As I mentioned in the OP, that is not the case for the covariance matrix V_0. We copy V_0 to V_0_new (line 321), update V_0_new (lines 365 and 397), but discard it all together.

I believe the funcion's documentation is already contradicted by the code and should not be considered as decisive evidence.

Nowcasting/functions/dfm.m

Lines 308 to 403 in 19f365c

	% Running the Kalman filter and smoother with current parameters
	% Note that log-liklihood is NOT re-estimated after the runKF step: This
	% effectively gives the previous iteration's log-likelihood
	% For more information on output, see runKF
	[Zsmooth, Vsmooth, VVsmooth, loglik] = runKF(y, A, C, Q, R, Z_0, V_0);
	%% MAXIMIZATION STEP (TRANSITION EQUATION)
	% See (Banbura & Modugno, 2010) for details.
	% Initialize output
	A_new = A;
	Q_new = Q;
	V_0_new = V_0;
	%%% 2A. UPDATE FACTOR PARAMETERS INDIVIDUALLY ----------------------------
	for i = 1:num_blocks % Loop for each block: factors are uncorrelated
	% SETUP INDEXING
	r_i = r(i); % r_i = 1 if block is loaded
	rp = r_i*p;
	rp1 = sum(r(1:i-1))*ppC;
	b_subset = rp1+1:rp1+rp; % Subset blocks: Helps for subsetting Zsmooth, Vsmooth
	t_start = rp1+1; % Transition matrix factor idx start
	t_end = rp1+r_i*ppC; % Transition matrix factor idx end
	% ESTIMATE FACTOR PORTION OF Q, A
	% Note: EZZ, EZZ_BB, EZZ_FB are parts of equations 6 and 8 in BM 2010
	% E[f_t*f_t' | Omega_T]
	EZZ = Zsmooth(b_subset, 2:end) * Zsmooth(b_subset, 2:end)'...
	+sum(Vsmooth(b_subset, b_subset, 2:end) ,3);
	% E[f_{t-1}*f_{t-1}' | Omega_T]
	EZZ_BB = Zsmooth(b_subset, 1:end-1)*Zsmooth(b_subset, 1:end-1)'...
	+sum(Vsmooth(b_subset, b_subset, 1:end-1), 3);
	% E[f_t*f_{t-1}' | Omega_T]
	EZZ_FB = Zsmooth(b_subset, 2:end)*Zsmooth(b_subset, 1:end-1)'...
	+sum(VVsmooth(b_subset, b_subset, :), 3);
	% Select transition matrix/covariance matrix for block i
	A_i = A(t_start:t_end, t_start:t_end);
	Q_i = Q(t_start:t_end, t_start:t_end);
	% Equation 6: Estimate VAR(p) for factor
	A_i(1:r_i,1:rp) = EZZ_FB(1:r_i,1:rp) * inv(EZZ_BB(1:rp,1:rp));
	% Equation 8: Covariance matrix of residuals of VAR
	Q_i(1:r_i,1:r_i) = (EZZ(1:r_i,1:r_i) - A_i(1:r_i,1:rp)* EZZ_FB(1:r_i,1:rp)') / T;
	% Place updated results in output matrix
	A_new(t_start:t_end, t_start:t_end) = A_i;
	Q_new(t_start:t_end, t_start:t_end) = Q_i;
	V_0_new(t_start:t_end, t_start:t_end) =...
	Vsmooth(t_start:t_end, t_start:t_end,1);
	end
	%%% 2B. UPDATING PARAMETERS FOR IDIOSYNCRATIC COMPONENT ------------------
	rp1 = sum(r)*ppC; % Col size of factor portion
	niM = sum(i_idio(1:nM)); % Number of monthly values
	t_start = rp1+1; % Start of idiosyncratic component index
	i_subset = t_start:rp1+niM; % Gives indices for monthly idiosyncratic component values
	% Below 3 estimate the idiosyncratic component (for eqns 6, 8 BM 2010)
	% E[f_t*f_t' | \Omega_T]
	EZZ = diag(diag(Zsmooth(t_start:end, 2:end) * Zsmooth(t_start:end, 2:end)'))...
	+ diag(diag(sum(Vsmooth(t_start:end, t_start:end, 2:end), 3)));
	% E[f_{t-1}*f_{t-1}' | \Omega_T]
	EZZ_BB = diag(diag(Zsmooth(t_start:end, 1:end-1)* Zsmooth(t_start:end, 1:end-1)'))...
	+ diag(diag(sum(Vsmooth(t_start:end, t_start:end, 1:end-1), 3)));
	% E[f_t*f_{t-1}' | \Omega_T]
	EZZ_FB = diag(diag(Zsmooth(t_start:end, 2:end)*Zsmooth(t_start:end, 1:end-1)'))...
	+ diag(diag(sum(VVsmooth(t_start:end, t_start:end, :), 3)));
	A_i = EZZ_FB * diag(1./diag((EZZ_BB))); % Equation 6
	Q_i = (EZZ - A_i*EZZ_FB') / T; % Equation 8
	% Place updated results in output matrix
	A_new(i_subset, i_subset) = A_i(1:niM,1:niM);
	Q_new(i_subset, i_subset) = Q_i(1:niM,1:niM);
	V_0_new(i_subset, i_subset) = diag(diag(Vsmooth(i_subset, i_subset, 1)));
	%% 3 MAXIMIZATION STEP (observation equation)
	%%% INITIALIZATION AND SETUP ----------------------------------------------
	Z_0 = Zsmooth(:,1); %zeros(size(Zsmooth,1),1); %

[MajesticKhan]

@victorbluhu No need to apologize!

Yes, I agree with you.
In fact I made a post about this earlier lol. Check this issue

@eric-qian @bbok can you please clarify this issue? Thank you!