Nothing
R/sparseDFM.R
In sparseDFM: Estimate Dynamic Factor Models with Sparse Loadings
Defines functions
sparseDFM
Documented in
sparseDFM
#' Estimate a Sparse Dynamic Factor Model
#'
#' @description
#' Main function to allow estimation of a DFM or a sparse DFM (with sparse loadings) on stationary data that may have arbitrary patterns of missing data. We allow the user:
#' \itemize{
#' \item an option for estimation method - \code{"PCA"}, \code{"2Stage"}, \code{"EM"} or \code{"EM-sparse"}
#' \item an option for \code{IID} or \code{AR1} idiosyncratic errors
#' \item an option for Kalman Filter/Smoother estimation using standard \code{multivariate} equations or fast \code{univariate} filtering equations
#' }
#'
#' @param X \code{n x p} numeric data matrix or data frame of (stationary) time series.
#' @param r Integer. Number of factors.
#' @param q Integer. The first q series (columns of X) should not be made sparse. Default q = 0.
#' @param alphas Numeric vector or value of LASSO regularisation parameters. Default is alphas = logspace(-2,3,100).
#' @param alg Character. Option for estimation algorithm. Default is \code{"EM-sparse"}. Options are:
#' \tabular{llll}{
#' \code{"PCA"} \tab\tab principle components analysis (PCA) for static factors seen in Stock and Watson (2002). \cr\cr
#' \code{"2Stage"} \tab\tab the two-stage framework of PCA plus Kalman filter/smoother seen in Giannone et al. (2008) and Doz et al. (2011). \cr\cr
#' \code{"EM"} \tab\tab the quasi-maximum likelihood approach using the EM algorithm to handle arbitrary patterns of missing data seen in Banbura and Modugno (2014). \cr\cr
#' \code{"EM-sparse"} \tab\tab the novel sparse EM approach allowing LASSO regularisation on factor loadings seen in (cite our paper). \cr\cr
#' }
#' @param err Character. Option for idiosyncratic errors. Default is \code{"IID"}. Options are:
#' \tabular{llll}{
#' \code{"IID"} \tab\tab errors are IID white noise. \cr\cr
#' \code{"AR1"} \tab\tab errors follow an AR(1) process. \cr\cr
#' }
#' @param kalman Character. Option for Kalman filter and smoother equations. Default is \code{"univariate"}. Options are:
#' \tabular{llll}{
#' \code{"multivariate"} \tab\tab classic Kalman filter and smoother equations seen in Shumway and Stoffer (1982). \cr\cr
#' \code{"univaraite"} \tab\tab univariate treatment (sequential processing) of the multivariate equations for fast Kalman filter and smoother seen in Koopman and Durbin (2000). \cr\cr
#' }
#' @param store.parameters Logical. Store outputs for every alpha L1 penalty parameter. Default is FALSE.
#' @param standardize Logical. Standardize the data before estimating the model. Default is \code{TRUE}.
#' @param max_iter Integer. Maximum number of EM iterations. Default is 100.
#' @param threshold Numeric. Tolerance on EM iterates. Default is 1e-4.
#'
#' @details
#' For full details of the model please refer to Mosley et al. (2023).
#'
#' @returns A list-of-lists-like S3 object of class 'sparseDFM' with the following elements:
#' \item{\code{data}}{A list containing information about the data with the following elements:
#' \tabular{llll}{
#' \code{X} \tab\tab is the original \eqn{n \times p}{n x p} numeric data matrix of (stationary) time series. \cr\cr
#' \code{standardize} \tab\tab is a logical value indicating whether the original data was standardized.\cr\cr
#' \code{X.mean} \tab\tab is a p-dimensional numeric vector of column means of \eqn{X}. \cr\cr
#' \code{X.sd} \tab\tab is a p-dimensional numeric vector of column standard deviations of \eqn{X}. \cr\cr
#' \code{X.bal} \tab\tab is a \eqn{n \times p}{n x p} numeric data matrix of the original \eqn{X} with missing data interpolated using \code{fillNA()}. \cr\cr
#' \code{eigen} \tab\tab is the eigen decomposition of \code{X.bal}. \cr\cr
#' \code{fitted} \tab\tab is the \eqn{n \times p}{n x p} predicted data matrix using the estimated parameters: \eqn{\hat{\Lambda}\hat{F}}{\hat{\Lambda}\hat{F}}. \cr\cr
#' \code{fitted.unscaled} \tab\tab is the \eqn{n \times p}{n x p} predicted data matrix using the estimated parameters: \eqn{\hat{\Lambda}\hat{F}}{\hat{\Lambda}\hat{F}} that has been unscaled back to original data scale if \code{standardize} is \code{TRUE}. \cr\cr
#' \code{method} \tab\tab the estimation algorithm used (\code{alg}). \cr\cr
#' \code{err} \tab\tab the type of idiosyncratic errors assumed. Either \code{IID} or \code{AR1}. \cr\cr
#' \code{call} \tab\tab call object obtained from \code{match.call()}. \cr\cr
#' }
#' }
#' \item{\code{params}}{A list containing the estimated parameters of the model with the following elements:
#' \tabular{llll}{
#' \code{A} \tab\tab the \eqn{r \times r}{r x r} factor transition matrix. \cr\cr
#' \code{Phi} \tab\tab the p-dimensional vector of AR(1) coefficients for the idiosyncratic errors. \cr\cr
#' \code{Lambda} \tab\tab the \eqn{p \times r}{p x r} loadings matrix. \cr\cr
#' \code{Sigma_u} \tab\tab the \eqn{r \times r}{r x r} factor transition error covariance matrix. \cr\cr
#' \code{Sigma_epsilon} \tab\tab the p-dimensional vector of idiosyncratic error variances. As \eqn{\bm{\Sigma}_{\epsilon}}{\Sigma(\epsilon)} is assumed to be diagonal. \cr\cr
#' }
#' }
#' \item{\code{state}}{A list containing the estimated states and state covariances with the following elements:
#' \tabular{llll}{
#' \code{factors} \tab\tab the \eqn{n \times r}{n x r} matrix of factor estimates. \cr\cr
#' \code{errors} \tab\tab the \eqn{n \times p}{n x p} matrix of AR(1) idiosyncratic error estimates. For err = AR1 only. \cr\cr
#' \code{factors.cov} \tab\tab the \eqn{r \times r \times n}{r x r x n} covariance matrices of the factor estimates. \cr\cr
#' \code{errors.cov} \tab\tab the \eqn{p \times p \times n}{p x p x n} covariance matrices of the AR(1) idiosyncratic error estimates. For err = AR1 only. \cr\cr
#' }
#' }
#' \item{\code{em}}{A list containing information about the EM algorithm with the following elements:
#' \tabular{llll}{
#' \code{converged} \tab\tab a logical value indicating whether the EM algorithm converged. \cr\cr
#' \code{alpha_grid} \tab\tab a numerical vector containing the LASSO tuning parameters considered in BIC evaluation before stopping. \cr\cr
#' \code{alpha_opt} \tab\tab the optimal LASSO tuning parameter used. \cr\cr
#' \code{bic} \tab\tab a numerical vector containing BIC values for the corresponding LASSO tuning parameter in \code{alpha_grid}. \cr\cr
#' \code{loglik} \tab\tab the log-likelihood of the innovations from the Kalman filter in the final model. \cr\cr
#' \code{num_iter} \tab\tab number of iterations taken by the EM algorithm. \cr\cr
#' \code{tol} \tab\tab tolerance for EM convergence. Matches \code{threshold} in the input. \cr\cr
#' \code{max_iter} \tab\tab maximum number of iterations allowed for the EM algorithm. Matches \code{max_iter} in the input. \cr\cr
#' \code{em_time} \tab\tab time taken for EM convergence \cr\cr
#' }
#' }
#' \item{\code{alpha.output}}{Parameter and state outputs for each L1-norm penalty parameter in \code{alphas} if \code{store.parameters = TRUE}.
#' }
#'
#'
#' @references
#' Banbura, M., & Modugno, M. (2014). Maximum likelihood estimation of factor models on datasets with arbitrary pattern of missing data. \emph{Journal of Applied Econometrics, 29}(1), 133-160.
#'
#' Doz, C., Giannone, D., & Reichlin, L. (2011). A two-step estimator for large approximate dynamic factor models based on Kalman filtering. \emph{Journal of Econometrics, 164}(1), 188-205.
#'
#' Giannone, D., Reichlin, L., & Small, D. (2008). Nowcasting: The real-time informational content of macroeconomic data. \emph{Journal of monetary economics, 55}(4), 665-676.
#'
#' Koopman, S. J., & Durbin, J. (2000). Fast filtering and smoothing for multivariate state space models. \emph{Journal of Time Series Analysis, 21}(3), 281-296.
#'
#' Mosley, L., Chan, TS., & Gibberd, A. (2023). sparseDFM: An R Package to Estimate Dynamic Factor Models with Sparse Loadings.
#'
#' Shumway, R. H., & Stoffer, D. S. (1982). An approach to time series smoothing and forecasting using the EM algorithm. \emph{Journal of time series analysis, 3}(4), 253-264.
#'
#' Stock, J. H., & Watson, M. W. (2002). Forecasting using principal components from a large number of predictors. \emph{Journal of the American statistical association, 97}(460), 1167-1179.
#'
#' @examples
#' # load inflation data set
#' data = inflation
#'
#' # reduce the size for these examples - full data found in vignette
#' data = data[1:60,]
#'
#' # make stationary by taking first differences
#' new_data = transformData(data, rep(2,ncol(data)))
#'
#' # tune for the number of factors to use
#' tuneFactors(new_data, type = 2)
#'
#' # fit a PCA using 3 PC's
#' fit.pca <- sparseDFM(new_data, r = 3, alg = 'PCA')
#'
#' # fit a DFM using the two-stage approach
#' fit.2stage <- sparseDFM(new_data, r = 3, alg = '2Stage')
#'
#' # fit a DFM using EM algorithm with 3 factors
#' fit.dfm <- sparseDFM(new_data, r = 3, alg = 'EM')
#'
#' # fit a Sparse DFM with 3 factors
#' fit.sdfm <- sparseDFM(new_data, r = 3, alg = 'EM-sparse')
#'
#' # observe the factor loadings of the sparse DFM
#' plot(fit.sdfm, type = 'loading.heatmap')
#'
#' # observe the factors
#' plot(fit.sdfm, type = 'factor')
#'
#' # observe the residuals
#' plot(fit.sdfm, type = 'residual')
#'
#' # observe the LASSO parameter selected and BIC values
#' plot(fit.sdfm, type = 'lasso.bic')
#'
#' # predict 3 steps ahead
#' predict(fit.sdfm, h = 3)
#'
#'
#' @useDynLib sparseDFM, .registration = TRUE
#'
#' @export
sparseDFM <- function(X, r, q = 0, alphas = logspace(-2,3,100), alg = 'EM-sparse', err = 'IID', kalman = 'univariate', store.parameters = FALSE, standardize = TRUE, max_iter=100, threshold=1e-4) {
## Correct input checks
if(alg != 'PCA' && alg != '2Stage' && alg != 'EM' && alg != 'EM-sparse'){
stop("Incorrect alg input")
}
if(err != 'IID' && err != 'AR1'){
stop("Incorrect err input")
}
if(kalman != 'multivariate' && kalman != 'univariate'){
stop("Incorrect kalman input")
}
if(!is.numeric(r) || r <= 0L){
stop("r needs to be an integer > 0")
}
if(!is.numeric(q) || q < 0){
stop("q needs to be an integer >= 0")
}
if(!is.numeric(alphas) || anyNA(alphas)){
stop("alphas must be a numeric vector with no missing values")
}
if(!is.numeric(max_iter) || max_iter <= 0){
stop("max_iter must be an integer > 0")
}
if(!is.numeric(threshold) || threshold <= 0){
stop("threshold must be > 0")
}
# return original X input
X.input = X
# Unclass X - make sure it is a numeric matrix
X = unclass(X)
# make sure lasso parameter grid is low to high
alphas = sort(alphas)
# dimensions
n = dim(X)[1]
p = dim(X)[2]
k = r + p
# standardize if TRUE
X.scale = scale(X)
X.mean = attr(X.scale, "scaled:center")
X.sd = attr(X.scale, "scaled:scale")
if(standardize){
X = X.scale
}
# Remove any columns that are entirely NA
numberNA = as.numeric(colSums(is.na(X)))
allNA = which(numberNA == n)
if(length(allNA) > 0){
X = X[,-allNA]
X.input = X.input[,-allNA]
p = dim(X)[2]
X.mean = X.mean[-allNA]
X.sd = X.sd[-allNA]
message('Columns: ',allNA, ' are entirely missing. Removing these columns.')
}
# label variables
obs.names = unlist(dimnames(X)[1])
series.names = unlist(dimnames(X)[2])
factor.names = paste0('F',1:r)
f.error.names = paste0('u',1:r)
## Apply algorithm: PCA, 2Stage, EM, EM-sparse
if(alg == 'PCA'){ # PCA algorithm applied
## PCA and VAR(1)
initialise <- initPCA(X,r,err)
a0_0 = initialise$a0_0
P0_0 = initialise$P0_0
A.tilde = initialise$A.tilde
Lambda.tilde = initialise$Lambda.tilde
Sigma.u.tilde = initialise$Sigma.u.tilde
Sigma.eta = initialise$Sigma.eta
factors.PCA = initialise$factors.pca
loadings.PCA = initialise$loadings.pca
fit_x = factors.PCA %*% t(loadings.PCA)
if(standardize){
fit_X = kronecker(t(X.sd),rep(1,n))*fit_x + kronecker(t(X.mean),rep(1,n))
}else{
fit_X = fit_x
}
dimnames(fit_x) = dimnames(X)
dimnames(fit_X) = dimnames(X)
dimnames(factors.PCA) = list(obs.names, factor.names)
## Output for PCA - depends on if err = 'AR1' or 'IID'
if(err == 'AR1'){
A = A.tilde[1:r,1:r]
Phi = A.tilde[(r+1):k,(r+1):k]
Lambda = Lambda.tilde[,1:r]
Sigma_u = Sigma.u.tilde[1:r,1:r]
Sigma_epsilon = Sigma.u.tilde[(r+1):k,(r+1):k]
factors.cov = P0_0[1:r,1:r]
dimnames(A) = list(factor.names,factor.names)
dimnames(Phi) = list(series.names,series.names)
dimnames(Lambda) = list(series.names, factor.names)
dimnames(Sigma_u) = list(f.error.names, f.error.names)
dimnames(Sigma_epsilon) = list(series.names, series.names)
dimnames(factors.cov) = list(factor.names, factor.names)
output = list(data = list(X = X.input,
standardize = standardize,
X.mean = X.mean,
X.sd = X.sd,
X.bal = initialise$X.bal,
eigen = initialise$eigen,
fitted = fit_x,
fitted.unscaled = fit_X,
method = alg,
err = err,
call = match.call()),
params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = A,
Phi = diag(Phi),
Lambda = Lambda,
Sigma_u = Sigma_u,
Sigma_epsilon = diag(Sigma_epsilon)),
state = list(factors = factors.PCA,
factors.cov = factors.cov))
}else {
A = A.tilde
Lambda = Lambda.tilde
Sigma_u = Sigma.u.tilde
Sigma_epsilon = Sigma.eta
factors.cov = P0_0[1:r,1:r]
dimnames(A) = list(factor.names,factor.names)
dimnames(Lambda) = list(series.names, factor.names)
dimnames(Sigma_u) = list(f.error.names, f.error.names)
dimnames(Sigma_epsilon) = list(series.names, series.names)
dimnames(factors.cov) = list(factor.names, factor.names)
output = list(data = list(X = X.input,
standardize = standardize,
X.mean = X.mean,
X.sd = X.sd,
X.bal = initialise$X.bal,
eigen = initialise$eigen,
fitted = fit_x,
fitted.unscaled = fit_X,
method = alg,
err = err,
call = match.call()),
params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = A,
Lambda = Lambda,
Sigma_u = Sigma_u,
Sigma_epsilon = diag(Sigma_epsilon)),
state = list(factors = factors.PCA,
factors.cov = factors.cov))
}
class(output) <- 'sparseDFM'
return(output)
}else if(alg == '2Stage'){ # 2 stage algorithm applied (Doz, 2011)
## Initialise with PCA and VAR(1)
initialise <- initPCA(X,r,err)
a0_0 = initialise$a0_0
P0_0 = initialise$P0_0
A.tilde = initialise$A.tilde
Lambda.tilde = initialise$Lambda.tilde
Sigma.u.tilde = initialise$Sigma.u.tilde
Sigma.eta = initialise$Sigma.eta
factors.PCA = initialise$factors.pca
loadings.PCA = initialise$loadings.pca
## Kalman Filter and Smoother (Doz (2011) 2-step)
if(kalman == 'univariate'){
KFS <- kalmanUnivariate(X, a0_0, P0_0, A.tilde, Lambda.tilde, Sigma.eta, Sigma.u.tilde)
}else{
KFS <- kalmanMultivariate(X, a0_0, P0_0, A.tilde, Lambda.tilde, Sigma.eta, Sigma.u.tilde)
}
state.KS = t(as.matrix(KFS$at_n))
covariance.KS = KFS$Pt_n
## Fill in missing data in X
fit_x = state.KS[,1:r] %*% t(Lambda.tilde[,1:r])
if(standardize){
fit_X = kronecker(t(X.sd),rep(1,n))*fit_x + kronecker(t(X.mean),rep(1,n))
}else{
fit_X = fit_x
}
dimnames(fit_x) = dimnames(X)
dimnames(fit_X) = dimnames(X)
## Output for 2Stage - depends on if err = 'AR1' or 'IID'
if(err == 'AR1'){
factors = state.KS[,1:r]
errors = state.KS[,(r+1):k]
dimnames(factors) = list(obs.names, factor.names)
dimnames(errors) = list(obs.names, series.names)
A = A.tilde[1:r,1:r]
Phi = A.tilde[(r+1):k,(r+1):k]
Lambda = Lambda.tilde[,1:r]
Sigma_u = Sigma.u.tilde[1:r,1:r]
Sigma_epsilon = Sigma.u.tilde[(r+1):k,(r+1):k]
factors.cov = covariance.KS[1:r,1:r,]
errors.cov = covariance.KS[(r+1):k,(r+1):k,]
dimnames(A) = list(factor.names,factor.names)
dimnames(Phi) = list(series.names,series.names)
dimnames(Lambda) = list(series.names, factor.names)
dimnames(Sigma_u) = list(f.error.names, f.error.names)
dimnames(Sigma_epsilon) = list(series.names, series.names)
dimnames(factors.cov) = list(factor.names, factor.names)
dimnames(errors.cov) = list(series.names, series.names)
output = list(data = list(X = X.input,
standardize = standardize,
X.mean = X.mean,
X.sd = X.sd,
X.bal = initialise$X.bal,
eigen = initialise$eigen,
fitted = fit_x,
fitted.unscaled = fit_X,
method = alg,
err = err,
call = match.call()),
params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = A,
Phi = diag(Phi),
Lambda = Lambda,
Sigma_u = Sigma_u,
Sigma_epsilon = diag(Sigma_epsilon)),
state = list(factors = factors,
errors = errors,
factors.cov = factors.cov,
errors.cov = errors.cov))
}else {
factors = state.KS
dimnames(factors) = list(obs.names, factor.names)
A = A.tilde
Lambda = Lambda.tilde
Sigma_u = Sigma.u.tilde
Sigma_epsilon = Sigma.eta
factors.cov = covariance.KS
dimnames(A) = list(factor.names,factor.names)
dimnames(Lambda) = list(series.names, factor.names)
dimnames(Sigma_u) = list(f.error.names, f.error.names)
dimnames(Sigma_epsilon) = list(series.names, series.names)
dimnames(factors.cov) = list(factor.names, factor.names)
output = list(data = list(X = X.input,
standardize = standardize,
X.mean = X.mean,
X.sd = X.sd,
X.bal = initialise$X.bal,
eigen = initialise$eigen,
fitted = fit_x,
fitted.unscaled = fit_X,
method = alg,
err = err,
call = match.call()),
params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = A,
Lambda = Lambda,
Sigma_u = Sigma_u,
Sigma_epsilon = diag(Sigma_epsilon)),
state = list(factors = factors,
factors.cov = factors.cov))
}
class(output) <- 'sparseDFM'
return(output)
}else if(alg == 'EM'){ # EM algorithm applied (Banbura and Modugno, 2014)
## Initialise with PCA and VAR(1)
initialise <- initPCA(X,r,err)
a0_0 = initialise$a0_0
P0_0 = initialise$P0_0
A.tilde = initialise$A.tilde
Lambda.tilde = initialise$Lambda.tilde
Sigma.u.tilde = initialise$Sigma.u.tilde
Sigma.eta = initialise$Sigma.eta
factors.PCA = initialise$factors.pca
loadings.PCA = initialise$loadings.pca
## EM Algorithm
# EM iterations function
st.em = Sys.time()
EM.fit <- EM(X, a0_0, P0_0, A.tilde, Lambda.tilde, Sigma.eta, Sigma.u.tilde,
err = err, kalman = kalman, sparse = FALSE, max_iter = max_iter, threshold = threshold)
et.em = Sys.time()
time.em = et.em - st.em
# Optimal parameters
a0_0 = EM.fit$a0_0
P0_0 = EM.fit$P0_0
A.tilde = EM.fit$A.tilde
Lambda.tilde = EM.fit$Lambda.tilde
Sigma.u.tilde = EM.fit$Sigma.u.tilde
Sigma.eta = EM.fit$Sigma.eta
loglik.store = EM.fit$loglik.store
converged = EM.fit$converged
num_iter = EM.fit$num_iter
# run KFS on final parameter estimates
if(kalman == 'univariate'){
best.KFS <- kalmanUnivariate(X, a0_0, P0_0, A.tilde, Lambda.tilde, Sigma.eta, Sigma.u.tilde)
}else{
best.KFS <- kalmanMultivariate(X, a0_0, P0_0, A.tilde, Lambda.tilde, Sigma.eta, Sigma.u.tilde)
}
state.EM = t(best.KFS$at_n)
covariance.EM = best.KFS$Pt_n
# fill in missing data in X
fit_x = state.EM[,1:r] %*% t(Lambda.tilde[,1:r])
if(standardize){
fit_X = kronecker(t(X.sd),rep(1,n))*fit_x + kronecker(t(X.mean),rep(1,n))
}else{
fit_X = fit_x
}
dimnames(fit_x) = dimnames(X)
dimnames(fit_X) = dimnames(X)
## Output for EM - depends on if err = 'AR1' or 'IID'
if(err == 'AR1'){
factors = state.EM[,1:r]
errors = state.EM[,(r+1):k]
dimnames(factors) = list(obs.names, factor.names)
dimnames(errors) = list(obs.names, series.names)
A = A.tilde[1:r,1:r]
Phi = A.tilde[(r+1):k,(r+1):k]
Lambda = Lambda.tilde[,1:r]
Sigma_u = Sigma.u.tilde[1:r,1:r]
Sigma_epsilon = Sigma.u.tilde[(r+1):k,(r+1):k]
factors.cov = covariance.EM[1:r,1:r,]
errors.cov = covariance.EM[(r+1):k,(r+1):k,]
dimnames(A) = list(factor.names,factor.names)
dimnames(Phi) = list(series.names,series.names)
dimnames(Lambda) = list(series.names, factor.names)
dimnames(Sigma_u) = list(f.error.names, f.error.names)
dimnames(Sigma_epsilon) = list(series.names, series.names)
dimnames(factors.cov) = list(factor.names, factor.names)
dimnames(errors.cov) = list(series.names, series.names)
output = list(data = list(X = X.input,
standardize = standardize,
X.mean = X.mean,
X.sd = X.sd,
X.bal = initialise$X.bal,
eigen = initialise$eigen,
fitted = fit_x,
fitted.unscaled = fit_X,
method = alg,
err = err,
call = match.call()),
params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = A,
Phi = diag(Phi),
Lambda = Lambda,
Sigma_u = Sigma_u,
Sigma_epsilon = diag(Sigma_epsilon)),
state = list(factors = factors,
errors = errors,
factors.cov = factors.cov,
errors.cov = errors.cov),
em = list(converged = converged,
loglik = loglik.store,
num_iter = num_iter,
tol = threshold,
max_iter = max_iter,
em_time = time.em))
}else {
factors = state.EM
dimnames(factors) = list(obs.names, factor.names)
A = A.tilde
Lambda = Lambda.tilde
Sigma_u = Sigma.u.tilde
Sigma_epsilon = Sigma.eta
factors.cov = covariance.EM
dimnames(A) = list(factor.names,factor.names)
dimnames(Lambda) = list(series.names, factor.names)
dimnames(Sigma_u) = list(f.error.names, f.error.names)
dimnames(Sigma_epsilon) = list(series.names, series.names)
dimnames(factors.cov) = list(factor.names, factor.names)
output = list(data = list(X = X.input,
standardize = standardize,
X.mean = X.mean,
X.sd = X.sd,
X.bal = initialise$X.bal,
eigen = initialise$eigen,
fitted = fit_x,
fitted.unscaled = fit_X,
method = alg,
err = err,
call = match.call()),
params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = A,
Lambda = Lambda,
Sigma_u = Sigma_u,
Sigma_epsilon = diag(Sigma_epsilon)),
state = list(factors = factors,
factors.cov = factors.cov),
em = list(converged = converged,
loglik = loglik.store,
num_iter = num_iter,
tol = threshold,
max_iter = max_iter,
em_time = time.em))
}
class(output) <- 'sparseDFM'
return(output)
}else { # sparse EM algorithm applied (Mosley et al, 2022)
## Initialise with PCA and VAR(1)
initialise <- initPCA(X,r,err)
a0_0 = initialise$a0_0
P0_0 = initialise$P0_0
A.tilde = initialise$A.tilde
Lambda.tilde = initialise$Lambda.tilde
Sigma.u.tilde = initialise$Sigma.u.tilde
Sigma.eta = initialise$Sigma.eta
factors.PCA = initialise$factors.pca
loadings.PCA = initialise$loadings.pca
## Sparsified EM Algorithm
## Loop over alphas and calculate BIC until column of Lambda becomes 0
bic <- c()
num_iter = c()
best.bic <- .Machine$double.xmax
alpha.output = list()
time.emsparse = c()
for(alphas.index in 1:length(alphas)) {
# lasso regularisation parameter
alpha.value = alphas[alphas.index]
# make into a matrix
alpha.value = matrix(alpha.value, nrow = p, ncol = r)
# alpha parameter set to 0 for no regularisation on first q series
if(q > 0){
alpha.value[1:q,] = 0
}
# EM iterations function
st.emsparse = Sys.time()
EM.fit <- EM(X, a0_0, P0_0, A.tilde, Lambda.tilde, Sigma.eta, Sigma.u.tilde,
alpha.lasso = alpha.value, err = err, kalman = kalman, sparse = TRUE,
max_iter = max_iter, threshold = threshold)
et.emsparse = Sys.time()
time.emsparse[alphas.index] = et.emsparse - st.emsparse
# store number of iterations for each alpha
num_iter[alphas.index] = EM.fit$num_iter
# check if a column of Lambda has been set entirely to 0
if(any(colSums(EM.fit$Lambda.tilde[(q+1):p,1:r]) == 0) && alphas.index == 1){
stop('First alpha value to large. All variables set to 0 in a factor. Make alpha value smaller.')
}
if(any(colSums(EM.fit$Lambda.tilde[(q+1):p,1:r]) == 0)){
break
}
# Update parameters - used for warm start of the EM algorithm
a0_0 = EM.fit$a0_0
P0_0 = EM.fit$P0_0
A.tilde = EM.fit$A.tilde
Lambda.tilde = EM.fit$Lambda.tilde
Sigma.u.tilde = EM.fit$Sigma.u.tilde
Sigma.eta = EM.fit$Sigma.eta
# run KFS on final parameter estimates
if(kalman == 'univariate'){
KFS <- kalmanUnivariate(X, a0_0, P0_0, A.tilde, Lambda.tilde, Sigma.eta, Sigma.u.tilde)
}else{
KFS <- kalmanMultivariate(X, a0_0, P0_0, A.tilde, Lambda.tilde, Sigma.eta, Sigma.u.tilde)
}
# calculate BIC
bic[alphas.index] = bicFunction(X, t(KFS$at_n[1:r,]), Lambda.tilde[,1:r])
# Store parameters for each alpha if store.parameters = TRUE
if(store.parameters){
store.state = t(KFS$at_n)
store.cov = KFS$Pt_n
if(err == 'AR1'){
store.factors = store.state[,1:r]
store.errors = store.state[,(r+1):k]
dimnames(store.factors) = list(obs.names, factor.names)
dimnames(store.errors) = list(obs.names, series.names)
store.A = A.tilde[1:r,1:r]
store.Phi = A.tilde[(r+1):k,(r+1):k]
store.Lambda = Lambda.tilde[,1:r]
store.Sigma_u = Sigma.u.tilde[1:r,1:r]
store.Sigma_epsilon = Sigma.u.tilde[(r+1):k,(r+1):k]
store.factors.cov = store.cov[1:r,1:r,]
store.errors.cov = store.cov[(r+1):k,(r+1):k,]
dimnames(store.A) = list(factor.names,factor.names)
dimnames(store.Phi) = list(series.names,series.names)
dimnames(store.Lambda) = list(series.names, factor.names)
dimnames(store.Sigma_u) = list(f.error.names, f.error.names)
dimnames(store.Sigma_epsilon) = list(series.names, series.names)
dimnames(store.factors.cov) = list(factor.names, factor.names)
dimnames(store.errors.cov) = list(series.names, series.names)
alpha.output[[alphas.index]] = list(params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = store.A,
Phi = diag(store.Phi),
Lambda = store.Lambda,
Sigma_u = store.Sigma_u,
Sigma_epsilon = diag(store.Sigma_epsilon)),
state = list(factors = store.factors,
errors = store.errors,
factors.cov = store.factors.cov,
errors.cov = store.errors.cov))
}else{
store.factors = store.state
dimnames(store.factors) = list(obs.names, factor.names)
store.A = A.tilde
store.Lambda = Lambda.tilde
store.Sigma_u = Sigma.u.tilde
store.Sigma_epsilon = Sigma.eta
store.factors.cov = store.cov
dimnames(store.A) = list(factor.names,factor.names)
dimnames(store.Lambda) = list(series.names, factor.names)
dimnames(store.Sigma_u) = list(f.error.names, f.error.names)
dimnames(store.Sigma_epsilon) = list(series.names, series.names)
dimnames(store.factors.cov) = list(factor.names, factor.names)
alpha.output[[alphas.index]] = list(params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = store.A,
Lambda = store.Lambda,
Sigma_u = store.Sigma_u,
Sigma_epsilon = diag(store.Sigma_epsilon)),
state = list(factors = store.factors,
factors.cov = store.factors.cov))
}
}
# store estimates if BIC improved
if(bic[alphas.index] < best.bic){
best.EM = EM.fit
best.KFS = KFS
best.bic = bic[alphas.index]
}
}
## Store the optimal outputs
time.emsparse = time.emsparse[1:length(bic)]
num_iter = num_iter[1:length(bic)]
alphas.used = alphas[1:length(bic)]
best.alpha = alphas[which.min(bic)]
loglik.store = best.EM$loglik.store
converged = best.EM$converged
a0_0 = best.EM$a0_0
P0_0 = best.EM$P0_0
A.tilde = best.EM$A.tilde
Lambda.tilde = best.EM$Lambda.tilde
Sigma.u.tilde = best.EM$Sigma.u.tilde
Sigma.eta = best.EM$Sigma.eta
state.EM = t(best.KFS$at_n)
covariance.EM = best.KFS$Pt_n
## Fill in missing data in X
fit_x = state.EM[,1:r] %*% t(Lambda.tilde[,1:r])
if(standardize){
fit_X = kronecker(t(X.sd),rep(1,n))*fit_x + kronecker(t(X.mean),rep(1,n))
}else{
fit_X = fit_x
}
dimnames(fit_x) = dimnames(X)
dimnames(fit_X) = dimnames(X)
if(store.parameters){
names(alpha.output) = paste0('alpha=',round(alphas.used,4))
}else{
alpha.output = NULL
}
## Output for EM-sparse - depends on if err = 'AR1' or 'IID'
if(err == 'AR1'){
factors = state.EM[,1:r]
errors = state.EM[,(r+1):k]
dimnames(factors) = list(obs.names, factor.names)
dimnames(errors) = list(obs.names, series.names)
A = A.tilde[1:r,1:r]
Phi = A.tilde[(r+1):k,(r+1):k]
Lambda = Lambda.tilde[,1:r]
Sigma_u = Sigma.u.tilde[1:r,1:r]
Sigma_epsilon = Sigma.u.tilde[(r+1):k,(r+1):k]
factors.cov = covariance.EM[1:r,1:r,]
errors.cov = covariance.EM[(r+1):k,(r+1):k,]
dimnames(A) = list(factor.names,factor.names)
dimnames(Phi) = list(series.names,series.names)
dimnames(Lambda) = list(series.names, factor.names)
dimnames(Sigma_u) = list(f.error.names, f.error.names)
dimnames(Sigma_epsilon) = list(series.names, series.names)
dimnames(factors.cov) = list(factor.names, factor.names)
dimnames(errors.cov) = list(series.names, series.names)
output = list(data = list(X = X.input,
standardize = standardize,
X.mean = X.mean,
X.sd = X.sd,
X.bal = initialise$X.bal,
eigen = initialise$eigen,
fitted = fit_x,
fitted.unscaled = fit_X,
method = alg,
err = err,
call = match.call()),
params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = A,
Phi = diag(Phi),
Lambda = Lambda,
Sigma_u = Sigma_u,
Sigma_epsilon = diag(Sigma_epsilon)),
state = list(factors = factors,
errors = errors,
factors.cov = factors.cov,
errors.cov = errors.cov),
em = list(converged = converged,
alpha_grid = alphas.used,
alpha_opt = best.alpha,
bic = bic,
loglik = loglik.store,
num_iter = num_iter,
tol = threshold,
max_iter = max_iter,
em_time = time.emsparse),
alpha.output = alpha.output)
}else {
factors = state.EM
dimnames(factors) = list(obs.names, factor.names)
A = A.tilde
Lambda = Lambda.tilde
Sigma_u = Sigma.u.tilde
Sigma_epsilon = Sigma.eta
factors.cov = covariance.EM
dimnames(A) = list(factor.names,factor.names)
dimnames(Lambda) = list(series.names, factor.names)
dimnames(Sigma_u) = list(f.error.names, f.error.names)
dimnames(Sigma_epsilon) = list(series.names, series.names)
dimnames(factors.cov) = list(factor.names, factor.names)
output = list(data = list(X = X.input,
standardize = standardize,
X.mean = X.mean,
X.sd = X.sd,
X.bal = initialise$X.bal,
eigen = initialise$eigen,
fitted = fit_x,
fitted.unscaled = fit_X,
method = alg,
err = err,
call = match.call()),
params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = A,
Lambda = Lambda,
Sigma_u = Sigma_u,
Sigma_epsilon = diag(Sigma_epsilon)),
state = list(factors = factors,
factors.cov = factors.cov),
em = list(converged = converged,
alpha_grid = alphas.used,
alpha_opt = best.alpha,
bic = bic,
loglik = loglik.store,
num_iter = num_iter,
tol = threshold,
max_iter = max_iter,
em_time = time.emsparse),
alpha.output = alpha.output)
}
class(output) <- 'sparseDFM'
return(output)
}
}
Try the sparseDFM package in your browser
Any scripts or data that you put into this service are public.
sparseDFM documentation built on March 31, 2023, 10:15 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions sparseDFM
Documented in sparseDFM
#' Estimate a Sparse Dynamic Factor Model
#'
#' @description
#' Main function to allow estimation of a DFM or a sparse DFM (with sparse loadings) on stationary data that may have arbitrary patterns of missing data. We allow the user:
#' \itemize{
#' \item an option for estimation method - \code{"PCA"}, \code{"2Stage"}, \code{"EM"} or \code{"EM-sparse"}
#' \item an option for \code{IID} or \code{AR1} idiosyncratic errors
#' \item an option for Kalman Filter/Smoother estimation using standard \code{multivariate} equations or fast \code{univariate} filtering equations
#' }
#'
#' @param X \code{n x p} numeric data matrix or data frame of (stationary) time series.
#' @param r Integer. Number of factors.
#' @param q Integer. The first q series (columns of X) should not be made sparse. Default q = 0.
#' @param alphas Numeric vector or value of LASSO regularisation parameters. Default is alphas = logspace(-2,3,100).
#' @param alg Character. Option for estimation algorithm. Default is \code{"EM-sparse"}. Options are:
#' \tabular{llll}{
#' \code{"PCA"} \tab\tab principle components analysis (PCA) for static factors seen in Stock and Watson (2002). \cr\cr
#' \code{"2Stage"} \tab\tab the two-stage framework of PCA plus Kalman filter/smoother seen in Giannone et al. (2008) and Doz et al. (2011). \cr\cr
#' \code{"EM"} \tab\tab the quasi-maximum likelihood approach using the EM algorithm to handle arbitrary patterns of missing data seen in Banbura and Modugno (2014). \cr\cr
#' \code{"EM-sparse"} \tab\tab the novel sparse EM approach allowing LASSO regularisation on factor loadings seen in (cite our paper). \cr\cr
#' }
#' @param err Character. Option for idiosyncratic errors. Default is \code{"IID"}. Options are:
#' \tabular{llll}{
#' \code{"IID"} \tab\tab errors are IID white noise. \cr\cr
#' \code{"AR1"} \tab\tab errors follow an AR(1) process. \cr\cr
#' }
#' @param kalman Character. Option for Kalman filter and smoother equations. Default is \code{"univariate"}. Options are:
#' \tabular{llll}{
#' \code{"multivariate"} \tab\tab classic Kalman filter and smoother equations seen in Shumway and Stoffer (1982). \cr\cr
#' \code{"univaraite"} \tab\tab univariate treatment (sequential processing) of the multivariate equations for fast Kalman filter and smoother seen in Koopman and Durbin (2000). \cr\cr
#' }
#' @param store.parameters Logical. Store outputs for every alpha L1 penalty parameter. Default is FALSE.
#' @param standardize Logical. Standardize the data before estimating the model. Default is \code{TRUE}.
#' @param max_iter Integer. Maximum number of EM iterations. Default is 100.
#' @param threshold Numeric. Tolerance on EM iterates. Default is 1e-4.
#'
#' @details
#' For full details of the model please refer to Mosley et al. (2023).
#'
#' @returns A list-of-lists-like S3 object of class 'sparseDFM' with the following elements:
#' \item{\code{data}}{A list containing information about the data with the following elements:
#' \tabular{llll}{
#' \code{X} \tab\tab is the original \eqn{n \times p}{n x p} numeric data matrix of (stationary) time series. \cr\cr
#' \code{standardize} \tab\tab is a logical value indicating whether the original data was standardized.\cr\cr
#' \code{X.mean} \tab\tab is a p-dimensional numeric vector of column means of \eqn{X}. \cr\cr
#' \code{X.sd} \tab\tab is a p-dimensional numeric vector of column standard deviations of \eqn{X}. \cr\cr
#' \code{X.bal} \tab\tab is a \eqn{n \times p}{n x p} numeric data matrix of the original \eqn{X} with missing data interpolated using \code{fillNA()}. \cr\cr
#' \code{eigen} \tab\tab is the eigen decomposition of \code{X.bal}. \cr\cr
#' \code{fitted} \tab\tab is the \eqn{n \times p}{n x p} predicted data matrix using the estimated parameters: \eqn{\hat{\Lambda}\hat{F}}{\hat{\Lambda}\hat{F}}. \cr\cr
#' \code{fitted.unscaled} \tab\tab is the \eqn{n \times p}{n x p} predicted data matrix using the estimated parameters: \eqn{\hat{\Lambda}\hat{F}}{\hat{\Lambda}\hat{F}} that has been unscaled back to original data scale if \code{standardize} is \code{TRUE}. \cr\cr
#' \code{method} \tab\tab the estimation algorithm used (\code{alg}). \cr\cr
#' \code{err} \tab\tab the type of idiosyncratic errors assumed. Either \code{IID} or \code{AR1}. \cr\cr
#' \code{call} \tab\tab call object obtained from \code{match.call()}. \cr\cr
#' }
#' }
#' \item{\code{params}}{A list containing the estimated parameters of the model with the following elements:
#' \tabular{llll}{
#' \code{A} \tab\tab the \eqn{r \times r}{r x r} factor transition matrix. \cr\cr
#' \code{Phi} \tab\tab the p-dimensional vector of AR(1) coefficients for the idiosyncratic errors. \cr\cr
#' \code{Lambda} \tab\tab the \eqn{p \times r}{p x r} loadings matrix. \cr\cr
#' \code{Sigma_u} \tab\tab the \eqn{r \times r}{r x r} factor transition error covariance matrix. \cr\cr
#' \code{Sigma_epsilon} \tab\tab the p-dimensional vector of idiosyncratic error variances. As \eqn{\bm{\Sigma}_{\epsilon}}{\Sigma(\epsilon)} is assumed to be diagonal. \cr\cr
#' }
#' }
#' \item{\code{state}}{A list containing the estimated states and state covariances with the following elements:
#' \tabular{llll}{
#' \code{factors} \tab\tab the \eqn{n \times r}{n x r} matrix of factor estimates. \cr\cr
#' \code{errors} \tab\tab the \eqn{n \times p}{n x p} matrix of AR(1) idiosyncratic error estimates. For err = AR1 only. \cr\cr
#' \code{factors.cov} \tab\tab the \eqn{r \times r \times n}{r x r x n} covariance matrices of the factor estimates. \cr\cr
#' \code{errors.cov} \tab\tab the \eqn{p \times p \times n}{p x p x n} covariance matrices of the AR(1) idiosyncratic error estimates. For err = AR1 only. \cr\cr
#' }
#' }
#' \item{\code{em}}{A list containing information about the EM algorithm with the following elements:
#' \tabular{llll}{
#' \code{converged} \tab\tab a logical value indicating whether the EM algorithm converged. \cr\cr
#' \code{alpha_grid} \tab\tab a numerical vector containing the LASSO tuning parameters considered in BIC evaluation before stopping. \cr\cr
#' \code{alpha_opt} \tab\tab the optimal LASSO tuning parameter used. \cr\cr
#' \code{bic} \tab\tab a numerical vector containing BIC values for the corresponding LASSO tuning parameter in \code{alpha_grid}. \cr\cr
#' \code{loglik} \tab\tab the log-likelihood of the innovations from the Kalman filter in the final model. \cr\cr
#' \code{num_iter} \tab\tab number of iterations taken by the EM algorithm. \cr\cr
#' \code{tol} \tab\tab tolerance for EM convergence. Matches \code{threshold} in the input. \cr\cr
#' \code{max_iter} \tab\tab maximum number of iterations allowed for the EM algorithm. Matches \code{max_iter} in the input. \cr\cr
#' \code{em_time} \tab\tab time taken for EM convergence \cr\cr
#' }
#' }
#' \item{\code{alpha.output}}{Parameter and state outputs for each L1-norm penalty parameter in \code{alphas} if \code{store.parameters = TRUE}.
#' }
#'
#'
#' @references
#' Banbura, M., & Modugno, M. (2014). Maximum likelihood estimation of factor models on datasets with arbitrary pattern of missing data. \emph{Journal of Applied Econometrics, 29}(1), 133-160.
#'
#' Doz, C., Giannone, D., & Reichlin, L. (2011). A two-step estimator for large approximate dynamic factor models based on Kalman filtering. \emph{Journal of Econometrics, 164}(1), 188-205.
#'
#' Giannone, D., Reichlin, L., & Small, D. (2008). Nowcasting: The real-time informational content of macroeconomic data. \emph{Journal of monetary economics, 55}(4), 665-676.
#'
#' Koopman, S. J., & Durbin, J. (2000). Fast filtering and smoothing for multivariate state space models. \emph{Journal of Time Series Analysis, 21}(3), 281-296.
#'
#' Mosley, L., Chan, TS., & Gibberd, A. (2023). sparseDFM: An R Package to Estimate Dynamic Factor Models with Sparse Loadings.
#'
#' Shumway, R. H., & Stoffer, D. S. (1982). An approach to time series smoothing and forecasting using the EM algorithm. \emph{Journal of time series analysis, 3}(4), 253-264.
#'
#' Stock, J. H., & Watson, M. W. (2002). Forecasting using principal components from a large number of predictors. \emph{Journal of the American statistical association, 97}(460), 1167-1179.
#'
#' @examples
#' # load inflation data set
#' data = inflation
#'
#' # reduce the size for these examples - full data found in vignette
#' data = data[1:60,]
#'
#' # make stationary by taking first differences
#' new_data = transformData(data, rep(2,ncol(data)))
#'
#' # tune for the number of factors to use
#' tuneFactors(new_data, type = 2)
#'
#' # fit a PCA using 3 PC's
#' fit.pca <- sparseDFM(new_data, r = 3, alg = 'PCA')
#'
#' # fit a DFM using the two-stage approach
#' fit.2stage <- sparseDFM(new_data, r = 3, alg = '2Stage')
#'
#' # fit a DFM using EM algorithm with 3 factors
#' fit.dfm <- sparseDFM(new_data, r = 3, alg = 'EM')
#'
#' # fit a Sparse DFM with 3 factors
#' fit.sdfm <- sparseDFM(new_data, r = 3, alg = 'EM-sparse')
#'
#' # observe the factor loadings of the sparse DFM
#' plot(fit.sdfm, type = 'loading.heatmap')
#'
#' # observe the factors
#' plot(fit.sdfm, type = 'factor')
#'
#' # observe the residuals
#' plot(fit.sdfm, type = 'residual')
#'
#' # observe the LASSO parameter selected and BIC values
#' plot(fit.sdfm, type = 'lasso.bic')
#'
#' # predict 3 steps ahead
#' predict(fit.sdfm, h = 3)
#'
#'
#' @useDynLib sparseDFM, .registration = TRUE
#'
#' @export
sparseDFM <- function(X, r, q = 0, alphas = logspace(-2,3,100), alg = 'EM-sparse', err = 'IID', kalman = 'univariate', store.parameters = FALSE, standardize = TRUE, max_iter=100, threshold=1e-4) {
## Correct input checks
if(alg != 'PCA' && alg != '2Stage' && alg != 'EM' && alg != 'EM-sparse'){
stop("Incorrect alg input")
}
if(err != 'IID' && err != 'AR1'){
stop("Incorrect err input")
}
if(kalman != 'multivariate' && kalman != 'univariate'){
stop("Incorrect kalman input")
}
if(!is.numeric(r) || r <= 0L){
stop("r needs to be an integer > 0")
}
if(!is.numeric(q) || q < 0){
stop("q needs to be an integer >= 0")
}
if(!is.numeric(alphas) || anyNA(alphas)){
stop("alphas must be a numeric vector with no missing values")
}
if(!is.numeric(max_iter) || max_iter <= 0){
stop("max_iter must be an integer > 0")
}
if(!is.numeric(threshold) || threshold <= 0){
stop("threshold must be > 0")
}
# return original X input
X.input = X
# Unclass X - make sure it is a numeric matrix
X = unclass(X)
# make sure lasso parameter grid is low to high
alphas = sort(alphas)
# dimensions
n = dim(X)[1]
p = dim(X)[2]
k = r + p
# standardize if TRUE
X.scale = scale(X)
X.mean = attr(X.scale, "scaled:center")
X.sd = attr(X.scale, "scaled:scale")
if(standardize){
X = X.scale
}
# Remove any columns that are entirely NA
numberNA = as.numeric(colSums(is.na(X)))
allNA = which(numberNA == n)
if(length(allNA) > 0){
X = X[,-allNA]
X.input = X.input[,-allNA]
p = dim(X)[2]
X.mean = X.mean[-allNA]
X.sd = X.sd[-allNA]
message('Columns: ',allNA, ' are entirely missing. Removing these columns.')
}
# label variables
obs.names = unlist(dimnames(X)[1])
series.names = unlist(dimnames(X)[2])
factor.names = paste0('F',1:r)
f.error.names = paste0('u',1:r)
## Apply algorithm: PCA, 2Stage, EM, EM-sparse
if(alg == 'PCA'){ # PCA algorithm applied
## PCA and VAR(1)
initialise <- initPCA(X,r,err)
a0_0 = initialise$a0_0
P0_0 = initialise$P0_0
A.tilde = initialise$A.tilde
Lambda.tilde = initialise$Lambda.tilde
Sigma.u.tilde = initialise$Sigma.u.tilde
Sigma.eta = initialise$Sigma.eta
factors.PCA = initialise$factors.pca
loadings.PCA = initialise$loadings.pca
fit_x = factors.PCA %*% t(loadings.PCA)
if(standardize){
fit_X = kronecker(t(X.sd),rep(1,n))*fit_x + kronecker(t(X.mean),rep(1,n))
}else{
fit_X = fit_x
}
dimnames(fit_x) = dimnames(X)
dimnames(fit_X) = dimnames(X)
dimnames(factors.PCA) = list(obs.names, factor.names)
## Output for PCA - depends on if err = 'AR1' or 'IID'
if(err == 'AR1'){
A = A.tilde[1:r,1:r]
Phi = A.tilde[(r+1):k,(r+1):k]
Lambda = Lambda.tilde[,1:r]
Sigma_u = Sigma.u.tilde[1:r,1:r]
Sigma_epsilon = Sigma.u.tilde[(r+1):k,(r+1):k]
factors.cov = P0_0[1:r,1:r]
dimnames(A) = list(factor.names,factor.names)
dimnames(Phi) = list(series.names,series.names)
dimnames(Lambda) = list(series.names, factor.names)
dimnames(Sigma_u) = list(f.error.names, f.error.names)
dimnames(Sigma_epsilon) = list(series.names, series.names)
dimnames(factors.cov) = list(factor.names, factor.names)
output = list(data = list(X = X.input,
standardize = standardize,
X.mean = X.mean,
X.sd = X.sd,
X.bal = initialise$X.bal,
eigen = initialise$eigen,
fitted = fit_x,
fitted.unscaled = fit_X,
method = alg,
err = err,
call = match.call()),
params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = A,
Phi = diag(Phi),
Lambda = Lambda,
Sigma_u = Sigma_u,
Sigma_epsilon = diag(Sigma_epsilon)),
state = list(factors = factors.PCA,
factors.cov = factors.cov))
}else {
A = A.tilde
Lambda = Lambda.tilde
Sigma_u = Sigma.u.tilde
Sigma_epsilon = Sigma.eta
factors.cov = P0_0[1:r,1:r]
dimnames(A) = list(factor.names,factor.names)
dimnames(Lambda) = list(series.names, factor.names)
dimnames(Sigma_u) = list(f.error.names, f.error.names)
dimnames(Sigma_epsilon) = list(series.names, series.names)
dimnames(factors.cov) = list(factor.names, factor.names)
output = list(data = list(X = X.input,
standardize = standardize,
X.mean = X.mean,
X.sd = X.sd,
X.bal = initialise$X.bal,
eigen = initialise$eigen,
fitted = fit_x,
fitted.unscaled = fit_X,
method = alg,
err = err,
call = match.call()),
params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = A,
Lambda = Lambda,
Sigma_u = Sigma_u,
Sigma_epsilon = diag(Sigma_epsilon)),
state = list(factors = factors.PCA,
factors.cov = factors.cov))
}
class(output) <- 'sparseDFM'
return(output)
}else if(alg == '2Stage'){ # 2 stage algorithm applied (Doz, 2011)
## Initialise with PCA and VAR(1)
initialise <- initPCA(X,r,err)
a0_0 = initialise$a0_0
P0_0 = initialise$P0_0
A.tilde = initialise$A.tilde
Lambda.tilde = initialise$Lambda.tilde
Sigma.u.tilde = initialise$Sigma.u.tilde
Sigma.eta = initialise$Sigma.eta
factors.PCA = initialise$factors.pca
loadings.PCA = initialise$loadings.pca
## Kalman Filter and Smoother (Doz (2011) 2-step)
if(kalman == 'univariate'){
KFS <- kalmanUnivariate(X, a0_0, P0_0, A.tilde, Lambda.tilde, Sigma.eta, Sigma.u.tilde)
}else{
KFS <- kalmanMultivariate(X, a0_0, P0_0, A.tilde, Lambda.tilde, Sigma.eta, Sigma.u.tilde)
}
state.KS = t(as.matrix(KFS$at_n))
covariance.KS = KFS$Pt_n
## Fill in missing data in X
fit_x = state.KS[,1:r] %*% t(Lambda.tilde[,1:r])
if(standardize){
fit_X = kronecker(t(X.sd),rep(1,n))*fit_x + kronecker(t(X.mean),rep(1,n))
}else{
fit_X = fit_x
}
dimnames(fit_x) = dimnames(X)
dimnames(fit_X) = dimnames(X)
## Output for 2Stage - depends on if err = 'AR1' or 'IID'
if(err == 'AR1'){
factors = state.KS[,1:r]
errors = state.KS[,(r+1):k]
dimnames(factors) = list(obs.names, factor.names)
dimnames(errors) = list(obs.names, series.names)
A = A.tilde[1:r,1:r]
Phi = A.tilde[(r+1):k,(r+1):k]
Lambda = Lambda.tilde[,1:r]
Sigma_u = Sigma.u.tilde[1:r,1:r]
Sigma_epsilon = Sigma.u.tilde[(r+1):k,(r+1):k]
factors.cov = covariance.KS[1:r,1:r,]
errors.cov = covariance.KS[(r+1):k,(r+1):k,]
dimnames(A) = list(factor.names,factor.names)
dimnames(Phi) = list(series.names,series.names)
dimnames(Lambda) = list(series.names, factor.names)
dimnames(Sigma_u) = list(f.error.names, f.error.names)
dimnames(Sigma_epsilon) = list(series.names, series.names)
dimnames(factors.cov) = list(factor.names, factor.names)
dimnames(errors.cov) = list(series.names, series.names)
output = list(data = list(X = X.input,
standardize = standardize,
X.mean = X.mean,
X.sd = X.sd,
X.bal = initialise$X.bal,
eigen = initialise$eigen,
fitted = fit_x,
fitted.unscaled = fit_X,
method = alg,
err = err,
call = match.call()),
params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = A,
Phi = diag(Phi),
Lambda = Lambda,
Sigma_u = Sigma_u,
Sigma_epsilon = diag(Sigma_epsilon)),
state = list(factors = factors,
errors = errors,
factors.cov = factors.cov,
errors.cov = errors.cov))
}else {
factors = state.KS
dimnames(factors) = list(obs.names, factor.names)
A = A.tilde
Lambda = Lambda.tilde
Sigma_u = Sigma.u.tilde
Sigma_epsilon = Sigma.eta
factors.cov = covariance.KS
dimnames(A) = list(factor.names,factor.names)
dimnames(Lambda) = list(series.names, factor.names)
dimnames(Sigma_u) = list(f.error.names, f.error.names)
dimnames(Sigma_epsilon) = list(series.names, series.names)
dimnames(factors.cov) = list(factor.names, factor.names)
output = list(data = list(X = X.input,
standardize = standardize,
X.mean = X.mean,
X.sd = X.sd,
X.bal = initialise$X.bal,
eigen = initialise$eigen,
fitted = fit_x,
fitted.unscaled = fit_X,
method = alg,
err = err,
call = match.call()),
params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = A,
Lambda = Lambda,
Sigma_u = Sigma_u,
Sigma_epsilon = diag(Sigma_epsilon)),
state = list(factors = factors,
factors.cov = factors.cov))
}
class(output) <- 'sparseDFM'
return(output)
}else if(alg == 'EM'){ # EM algorithm applied (Banbura and Modugno, 2014)
## Initialise with PCA and VAR(1)
initialise <- initPCA(X,r,err)
a0_0 = initialise$a0_0
P0_0 = initialise$P0_0
A.tilde = initialise$A.tilde
Lambda.tilde = initialise$Lambda.tilde
Sigma.u.tilde = initialise$Sigma.u.tilde
Sigma.eta = initialise$Sigma.eta
factors.PCA = initialise$factors.pca
loadings.PCA = initialise$loadings.pca
## EM Algorithm
# EM iterations function
st.em = Sys.time()
EM.fit <- EM(X, a0_0, P0_0, A.tilde, Lambda.tilde, Sigma.eta, Sigma.u.tilde,
err = err, kalman = kalman, sparse = FALSE, max_iter = max_iter, threshold = threshold)
et.em = Sys.time()
time.em = et.em - st.em
# Optimal parameters
a0_0 = EM.fit$a0_0
P0_0 = EM.fit$P0_0
A.tilde = EM.fit$A.tilde
Lambda.tilde = EM.fit$Lambda.tilde
Sigma.u.tilde = EM.fit$Sigma.u.tilde
Sigma.eta = EM.fit$Sigma.eta
loglik.store = EM.fit$loglik.store
converged = EM.fit$converged
num_iter = EM.fit$num_iter
# run KFS on final parameter estimates
if(kalman == 'univariate'){
best.KFS <- kalmanUnivariate(X, a0_0, P0_0, A.tilde, Lambda.tilde, Sigma.eta, Sigma.u.tilde)
}else{
best.KFS <- kalmanMultivariate(X, a0_0, P0_0, A.tilde, Lambda.tilde, Sigma.eta, Sigma.u.tilde)
}
state.EM = t(best.KFS$at_n)
covariance.EM = best.KFS$Pt_n
# fill in missing data in X
fit_x = state.EM[,1:r] %*% t(Lambda.tilde[,1:r])
if(standardize){
fit_X = kronecker(t(X.sd),rep(1,n))*fit_x + kronecker(t(X.mean),rep(1,n))
}else{
fit_X = fit_x
}
dimnames(fit_x) = dimnames(X)
dimnames(fit_X) = dimnames(X)
## Output for EM - depends on if err = 'AR1' or 'IID'
if(err == 'AR1'){
factors = state.EM[,1:r]
errors = state.EM[,(r+1):k]
dimnames(factors) = list(obs.names, factor.names)
dimnames(errors) = list(obs.names, series.names)
A = A.tilde[1:r,1:r]
Phi = A.tilde[(r+1):k,(r+1):k]
Lambda = Lambda.tilde[,1:r]
Sigma_u = Sigma.u.tilde[1:r,1:r]
Sigma_epsilon = Sigma.u.tilde[(r+1):k,(r+1):k]
factors.cov = covariance.EM[1:r,1:r,]
errors.cov = covariance.EM[(r+1):k,(r+1):k,]
dimnames(A) = list(factor.names,factor.names)
dimnames(Phi) = list(series.names,series.names)
dimnames(Lambda) = list(series.names, factor.names)
dimnames(Sigma_u) = list(f.error.names, f.error.names)
dimnames(Sigma_epsilon) = list(series.names, series.names)
dimnames(factors.cov) = list(factor.names, factor.names)
dimnames(errors.cov) = list(series.names, series.names)
output = list(data = list(X = X.input,
standardize = standardize,
X.mean = X.mean,
X.sd = X.sd,
X.bal = initialise$X.bal,
eigen = initialise$eigen,
fitted = fit_x,
fitted.unscaled = fit_X,
method = alg,
err = err,
call = match.call()),
params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = A,
Phi = diag(Phi),
Lambda = Lambda,
Sigma_u = Sigma_u,
Sigma_epsilon = diag(Sigma_epsilon)),
state = list(factors = factors,
errors = errors,
factors.cov = factors.cov,
errors.cov = errors.cov),
em = list(converged = converged,
loglik = loglik.store,
num_iter = num_iter,
tol = threshold,
max_iter = max_iter,
em_time = time.em))
}else {
factors = state.EM
dimnames(factors) = list(obs.names, factor.names)
A = A.tilde
Lambda = Lambda.tilde
Sigma_u = Sigma.u.tilde
Sigma_epsilon = Sigma.eta
factors.cov = covariance.EM
dimnames(A) = list(factor.names,factor.names)
dimnames(Lambda) = list(series.names, factor.names)
dimnames(Sigma_u) = list(f.error.names, f.error.names)
dimnames(Sigma_epsilon) = list(series.names, series.names)
dimnames(factors.cov) = list(factor.names, factor.names)
output = list(data = list(X = X.input,
standardize = standardize,
X.mean = X.mean,
X.sd = X.sd,
X.bal = initialise$X.bal,
eigen = initialise$eigen,
fitted = fit_x,
fitted.unscaled = fit_X,
method = alg,
err = err,
call = match.call()),
params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = A,
Lambda = Lambda,
Sigma_u = Sigma_u,
Sigma_epsilon = diag(Sigma_epsilon)),
state = list(factors = factors,
factors.cov = factors.cov),
em = list(converged = converged,
loglik = loglik.store,
num_iter = num_iter,
tol = threshold,
max_iter = max_iter,
em_time = time.em))
}
class(output) <- 'sparseDFM'
return(output)
}else { # sparse EM algorithm applied (Mosley et al, 2022)
## Initialise with PCA and VAR(1)
initialise <- initPCA(X,r,err)
a0_0 = initialise$a0_0
P0_0 = initialise$P0_0
A.tilde = initialise$A.tilde
Lambda.tilde = initialise$Lambda.tilde
Sigma.u.tilde = initialise$Sigma.u.tilde
Sigma.eta = initialise$Sigma.eta
factors.PCA = initialise$factors.pca
loadings.PCA = initialise$loadings.pca
## Sparsified EM Algorithm
## Loop over alphas and calculate BIC until column of Lambda becomes 0
bic <- c()
num_iter = c()
best.bic <- .Machine$double.xmax
alpha.output = list()
time.emsparse = c()
for(alphas.index in 1:length(alphas)) {
# lasso regularisation parameter
alpha.value = alphas[alphas.index]
# make into a matrix
alpha.value = matrix(alpha.value, nrow = p, ncol = r)
# alpha parameter set to 0 for no regularisation on first q series
if(q > 0){
alpha.value[1:q,] = 0
}
# EM iterations function
st.emsparse = Sys.time()
EM.fit <- EM(X, a0_0, P0_0, A.tilde, Lambda.tilde, Sigma.eta, Sigma.u.tilde,
alpha.lasso = alpha.value, err = err, kalman = kalman, sparse = TRUE,
max_iter = max_iter, threshold = threshold)
et.emsparse = Sys.time()
time.emsparse[alphas.index] = et.emsparse - st.emsparse
# store number of iterations for each alpha
num_iter[alphas.index] = EM.fit$num_iter
# check if a column of Lambda has been set entirely to 0
if(any(colSums(EM.fit$Lambda.tilde[(q+1):p,1:r]) == 0) && alphas.index == 1){
stop('First alpha value to large. All variables set to 0 in a factor. Make alpha value smaller.')
}
if(any(colSums(EM.fit$Lambda.tilde[(q+1):p,1:r]) == 0)){
break
}
# Update parameters - used for warm start of the EM algorithm
a0_0 = EM.fit$a0_0
P0_0 = EM.fit$P0_0
A.tilde = EM.fit$A.tilde
Lambda.tilde = EM.fit$Lambda.tilde
Sigma.u.tilde = EM.fit$Sigma.u.tilde
Sigma.eta = EM.fit$Sigma.eta
# run KFS on final parameter estimates
if(kalman == 'univariate'){
KFS <- kalmanUnivariate(X, a0_0, P0_0, A.tilde, Lambda.tilde, Sigma.eta, Sigma.u.tilde)
}else{
KFS <- kalmanMultivariate(X, a0_0, P0_0, A.tilde, Lambda.tilde, Sigma.eta, Sigma.u.tilde)
}
# calculate BIC
bic[alphas.index] = bicFunction(X, t(KFS$at_n[1:r,]), Lambda.tilde[,1:r])
# Store parameters for each alpha if store.parameters = TRUE
if(store.parameters){
store.state = t(KFS$at_n)
store.cov = KFS$Pt_n
if(err == 'AR1'){
store.factors = store.state[,1:r]
store.errors = store.state[,(r+1):k]
dimnames(store.factors) = list(obs.names, factor.names)
dimnames(store.errors) = list(obs.names, series.names)
store.A = A.tilde[1:r,1:r]
store.Phi = A.tilde[(r+1):k,(r+1):k]
store.Lambda = Lambda.tilde[,1:r]
store.Sigma_u = Sigma.u.tilde[1:r,1:r]
store.Sigma_epsilon = Sigma.u.tilde[(r+1):k,(r+1):k]
store.factors.cov = store.cov[1:r,1:r,]
store.errors.cov = store.cov[(r+1):k,(r+1):k,]
dimnames(store.A) = list(factor.names,factor.names)
dimnames(store.Phi) = list(series.names,series.names)
dimnames(store.Lambda) = list(series.names, factor.names)
dimnames(store.Sigma_u) = list(f.error.names, f.error.names)
dimnames(store.Sigma_epsilon) = list(series.names, series.names)
dimnames(store.factors.cov) = list(factor.names, factor.names)
dimnames(store.errors.cov) = list(series.names, series.names)
alpha.output[[alphas.index]] = list(params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = store.A,
Phi = diag(store.Phi),
Lambda = store.Lambda,
Sigma_u = store.Sigma_u,
Sigma_epsilon = diag(store.Sigma_epsilon)),
state = list(factors = store.factors,
errors = store.errors,
factors.cov = store.factors.cov,
errors.cov = store.errors.cov))
}else{
store.factors = store.state
dimnames(store.factors) = list(obs.names, factor.names)
store.A = A.tilde
store.Lambda = Lambda.tilde
store.Sigma_u = Sigma.u.tilde
store.Sigma_epsilon = Sigma.eta
store.factors.cov = store.cov
dimnames(store.A) = list(factor.names,factor.names)
dimnames(store.Lambda) = list(series.names, factor.names)
dimnames(store.Sigma_u) = list(f.error.names, f.error.names)
dimnames(store.Sigma_epsilon) = list(series.names, series.names)
dimnames(store.factors.cov) = list(factor.names, factor.names)
alpha.output[[alphas.index]] = list(params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = store.A,
Lambda = store.Lambda,
Sigma_u = store.Sigma_u,
Sigma_epsilon = diag(store.Sigma_epsilon)),
state = list(factors = store.factors,
factors.cov = store.factors.cov))
}
}
# store estimates if BIC improved
if(bic[alphas.index] < best.bic){
best.EM = EM.fit
best.KFS = KFS
best.bic = bic[alphas.index]
}
}
## Store the optimal outputs
time.emsparse = time.emsparse[1:length(bic)]
num_iter = num_iter[1:length(bic)]
alphas.used = alphas[1:length(bic)]
best.alpha = alphas[which.min(bic)]
loglik.store = best.EM$loglik.store
converged = best.EM$converged
a0_0 = best.EM$a0_0
P0_0 = best.EM$P0_0
A.tilde = best.EM$A.tilde
Lambda.tilde = best.EM$Lambda.tilde
Sigma.u.tilde = best.EM$Sigma.u.tilde
Sigma.eta = best.EM$Sigma.eta
state.EM = t(best.KFS$at_n)
covariance.EM = best.KFS$Pt_n
## Fill in missing data in X
fit_x = state.EM[,1:r] %*% t(Lambda.tilde[,1:r])
if(standardize){
fit_X = kronecker(t(X.sd),rep(1,n))*fit_x + kronecker(t(X.mean),rep(1,n))
}else{
fit_X = fit_x
}
dimnames(fit_x) = dimnames(X)
dimnames(fit_X) = dimnames(X)
if(store.parameters){
names(alpha.output) = paste0('alpha=',round(alphas.used,4))
}else{
alpha.output = NULL
}
## Output for EM-sparse - depends on if err = 'AR1' or 'IID'
if(err == 'AR1'){
factors = state.EM[,1:r]
errors = state.EM[,(r+1):k]
dimnames(factors) = list(obs.names, factor.names)
dimnames(errors) = list(obs.names, series.names)
A = A.tilde[1:r,1:r]
Phi = A.tilde[(r+1):k,(r+1):k]
Lambda = Lambda.tilde[,1:r]
Sigma_u = Sigma.u.tilde[1:r,1:r]
Sigma_epsilon = Sigma.u.tilde[(r+1):k,(r+1):k]
factors.cov = covariance.EM[1:r,1:r,]
errors.cov = covariance.EM[(r+1):k,(r+1):k,]
dimnames(A) = list(factor.names,factor.names)
dimnames(Phi) = list(series.names,series.names)
dimnames(Lambda) = list(series.names, factor.names)
dimnames(Sigma_u) = list(f.error.names, f.error.names)
dimnames(Sigma_epsilon) = list(series.names, series.names)
dimnames(factors.cov) = list(factor.names, factor.names)
dimnames(errors.cov) = list(series.names, series.names)
output = list(data = list(X = X.input,
standardize = standardize,
X.mean = X.mean,
X.sd = X.sd,
X.bal = initialise$X.bal,
eigen = initialise$eigen,
fitted = fit_x,
fitted.unscaled = fit_X,
method = alg,
err = err,
call = match.call()),
params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = A,
Phi = diag(Phi),
Lambda = Lambda,
Sigma_u = Sigma_u,
Sigma_epsilon = diag(Sigma_epsilon)),
state = list(factors = factors,
errors = errors,
factors.cov = factors.cov,
errors.cov = errors.cov),
em = list(converged = converged,
alpha_grid = alphas.used,
alpha_opt = best.alpha,
bic = bic,
loglik = loglik.store,
num_iter = num_iter,
tol = threshold,
max_iter = max_iter,
em_time = time.emsparse),
alpha.output = alpha.output)
}else {
factors = state.EM
dimnames(factors) = list(obs.names, factor.names)
A = A.tilde
Lambda = Lambda.tilde
Sigma_u = Sigma.u.tilde
Sigma_epsilon = Sigma.eta
factors.cov = covariance.EM
dimnames(A) = list(factor.names,factor.names)
dimnames(Lambda) = list(series.names, factor.names)
dimnames(Sigma_u) = list(f.error.names, f.error.names)
dimnames(Sigma_epsilon) = list(series.names, series.names)
dimnames(factors.cov) = list(factor.names, factor.names)
output = list(data = list(X = X.input,
standardize = standardize,
X.mean = X.mean,
X.sd = X.sd,
X.bal = initialise$X.bal,
eigen = initialise$eigen,
fitted = fit_x,
fitted.unscaled = fit_X,
method = alg,
err = err,
call = match.call()),
params = list(a0_0 = a0_0,
P0_0 = P0_0,
A = A,
Lambda = Lambda,
Sigma_u = Sigma_u,
Sigma_epsilon = diag(Sigma_epsilon)),
state = list(factors = factors,
factors.cov = factors.cov),
em = list(converged = converged,
alpha_grid = alphas.used,
alpha_opt = best.alpha,
bic = bic,
loglik = loglik.store,
num_iter = num_iter,
tol = threshold,
max_iter = max_iter,
em_time = time.emsparse),
alpha.output = alpha.output)
}
class(output) <- 'sparseDFM'
return(output)
}
}
Try the sparseDFM package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.