dfosr: MCMC Sampling Algorithm for the Dynamic Function-on-Scalars...
In drkowal/dfosr: Dynamic Function-on-Scalars Regression

Description Usage Arguments Value Note Examples

Runs the MCMC for the dynamic function-on-scalars regression model based on an reduced-rank expansion. The dynamic regression coefficients are modeled as random walks, with AR(1) models for the errors. Various shrinkage priors are available for the dynamic coefficient innovation terms.

dfosr(Y, tau, X = NULL, K = NULL, nsave = 1000, nburn = 1000,
  nskip = 3, mcmc_params = list("beta", "fk", "alpha", "mu_k",
  "ar_phi"), X_Tp1 = 1, use_obs_SV = FALSE,
  includeBasisInnovation = FALSE, Con_mat = NULL)

`Y`	the `T x m` data observation matrix, where `T` is the number of time points and `m` is the number of observation points (`NA`s allowed)
`tau`	the `m x d` matrix of coordinates of observation points
`X`	the `T x p` matrix of predictors; if NULL, only include an intercept
`K`	the number of factors; if NULL, use SVD-based proportion of variability explained
`nsave`	number of MCMC iterations to record
`nburn`	number of MCMC iterations to discard (burin-in)
`nskip`	number of MCMC iterations to skip between saving iterations, i.e., save every (nskip + 1)th draw
`mcmc_params`	named list of parameters for which we store the MCMC output; must be one or more of "beta" (dynamic factors) "fk" (loading curves) "alpha" (regression coefficients; possibly dynamic) "mu_k" (intercept term for factor k) "ar_phi" (AR coefficients for each k under AR(1) model) "sigma_et" (observation error SD; possibly dynamic) "Yhat" (fitted values) "Ypred" (posterior predictive values) "Yfore" (one-step forecast; includes the estimate and the distribution)
`X_Tp1`	the `p x 1` matrix of predictors at the forecasting time point `T + 1`
`use_obs_SV`	logical; when TRUE, include a stochastic volatility model for the observation error variance
`includeBasisInnovation`	logical; when TRUE, include an iid basis coefficient term for residual correlation (i.e., the idiosyncratic error term for a factor model on the full basis matrix)
`Con_mat`	a `Jc x m` matrix of constraints for the loading curves such that `Con_mat*fk = 0` for each loading curve `fk`; default is NULL for no constraints.

A named list of the nsave MCMC samples for the parameters named in mcmc_params

This sampler loops over the k=1,...,K factors, so the sampler is O(T*K*p^3) instead of O(T*(K*p)^3).

## Not run: 
# Simulate some data:
sim_data = simulate_dfosr(T = 200, m = 50, p_0 = 2, p_1 = 2)
Y = sim_data$Y; X = sim_data$X; tau = sim_data$tau
T = nrow(Y); m = ncol(Y); p = ncol(X) # Dimensions

# Run the MCMC w/ K = 6:
out = dfosr(Y = Y, tau = tau, X = X, K = 6,
           mcmc_params = list("beta", "fk", "alpha", "Yhat", "Ypred"))

# Plot a dynamic regression coefficient function
j = 3 # choose a predictor
post_alpha_tilde_j = get_post_alpha_tilde(out$fk, out$alpha[,,j,])

# Posterior mean:
alpha_tilde_j_pm = colMeans(post_alpha_tilde_j)

# Lower and Upper 95% credible intervals:
alpha_tilde_j_lower = apply(post_alpha_tilde_j, 2:3, quantile, c(0.05/2))
alpha_tilde_j_upper = apply(post_alpha_tilde_j, 2:3, quantile, c(1 - 0.05/2))

# Plot lower pointwise interval:
filled.contour(1:T, tau, alpha_tilde_j_lower,
               zlim = range(alpha_tilde_j_lower, alpha_tilde_j_upper),
               color = terrain.colors,
               xlab = 'Time', ylab = expression(tau),
               main = paste('Lower 95% Credible Intervals, j =',j))
# Plot posterior Mean:
filled.contour(1:T, tau, alpha_tilde_j_pm,
               zlim = range(alpha_tilde_j_lower, alpha_tilde_j_upper),
               color = terrain.colors,
               xlab = 'Time', ylab = expression(tau),
               main = paste('Posterior Mean, j =',j))
# Plot upper pointwise interval:
filled.contour(1:T, tau, sim_data$alpha_tilde_true[,j,],
               zlim = range(alpha_tilde_j_lower, alpha_tilde_j_upper),
               color = terrain.colors,
               xlab = 'Time', ylab = expression(tau),
               main = paste('Upper 95% Credible Intervals, j =',j))
# Truth:
filled.contour(1:T, tau, alpha_tilde_j_upper,
               zlim = range(alpha_tilde_j_lower, alpha_tilde_j_upper),
               color = terrain.colors,
               xlab = 'Time', ylab = expression(tau),
               main = paste('True regression coefficients, j =',j))

# Verify by plotting at two time slices:
t1 = ceiling(0.2*T); # Time t1
plot_curve(post_f = post_alpha_tilde_j[,t1,],
           tau = tau,
           main = paste('Predictor j =',j,'at time t =',t1))
# Add the true regression coefficient function:
lines(tau, sim_data$alpha_tilde_true[t1,j,], lwd=8, col='black', lty=6)

t2 = ceiling(0.8*T) # Time t2
plot_curve(post_f = post_alpha_tilde_j[,t2,],
           tau = tau,
           main = paste('Predictor j =',j,'at time t =',t2))
# Add the true regression coefficient function:
lines(tau, sim_data$alpha_tilde_true[t2,j,], lwd=8, col='black', lty=6)

# Plot the factors:
plot_factors(post_beta = out$beta)

# Plot the loading curves:
plot_flc(post_fk = out$fk, tau = tau)

# Plot a fitted value w/ posterior predictive credible intervals:
i = sample(1:T, 1); # Select a random time i
plot_fitted(y = Y[i,],
            mu = colMeans(out$Yhat)[i,],
            postY = out$Ypred[,i,],
            y_true = sim_data$Y_true[i,],
            t01 = tau)

## End(Not run)