fix predict for svm

NEWS

@@ -1,3 +1,9 @@
+bssm 1.1.4 (Release date: 2021-03-)
+==============
+   * Better documentation for SV model, and changed ordering of arguments to emphasise the 
+     recommended parameterization.
+   * Fixed predict method for SV model.
+
 bssm 1.1.3-2 (Release date: 2021-02-24)
 ==============
    * Fixed missing parenthesis causing compilation fail in case of no OpenMP support.

R/models.R

@@ -1258,12 +1258,13 @@ bsm_ng <- function(y, sd_level, sd_slope, sd_seasonal, sd_noise,
 #' first order autoregressive signal.
 #'
 #' @param y Vector or a \code{\link{ts}} object of observations.
-#' @param rho prior for autoregressive coefficient.
-#' @param sigma Prior for sigma parameter of observation equation.
 #' @param mu Prior for mu parameter of transition equation.
-#' Ignored if \code{sigma} is provided.
+#' @param rho prior for autoregressive coefficient.
 #' @param sd_ar Prior for the standard deviation of noise of the AR-process.
-#' @return Object of class \code{svm} or \code{svm2}.
+#' @param sigma Prior for sigma parameter of observation equation, internally denoted as phi. Ignored 
+#' if \code{mu} is provided. Note that typically parametrization using mu is preferred due to 
+#' better numerical properties and availability of better Gaussian approximation.
+#' @return Object of class \code{svm}.
 #' @export
 #' @rdname svm
 #' @examples
@@ -1284,7 +1285,7 @@ bsm_ng <- function(y, sd_level, sd_slope, sd_seasonal, sd_noise,
 #'   sd_ar = halfnormal(pars[2],sd=5),
 #'   sigma = halfnormal(pars[3],sd=2))
 #'
-svm <- function(y, rho, sd_ar, sigma, mu) {
+svm <- function(y, mu, rho, sd_ar, sigma) {
 
   if(!missing(sigma) && !missing(mu)) {
     stop("Define either sigma or mu, but not both.")
@@ -1294,7 +1295,7 @@ svm <- function(y, rho, sd_ar, sigma, mu) {
 
   check_rho(rho$init)
   check_sd(sd_ar$init, "rho")
-  if(missing(sigma)) {
+  if(!missing(mu)) {
     svm_type <- 1L
     check_mu(mu$init)
     initial_mode <- matrix(log(pmax(1e-4, y^2)), ncol = 1)

R/run_mcmc.R

@@ -258,6 +258,14 @@ run_mcmc.gaussian <- function(model, iter, output_type = "full",
 #'   theme(legend.title = element_blank()) + 
 #'   xlab("Time") + ylab("Level")
 #' 
+#' # theta
+#' d_theta <- as.data.frame(fit, variable = "theta")
+#' ggplot(d_theta, aes(x = value)) + 
+#'  geom_density(aes(weight = weight), adjust = 2, fill = "#92f0a8") + 
+#'  facet_wrap(~ variable, scales = "free") + 
+#'  theme_bw()
+#'  
+#'  
 #' # Bivariate Poisson model:
 #' 
 #' set.seed(1)

src/model_ssm_ung.cpp

@@ -392,7 +392,7 @@ arma::vec ssm_ung::log_weights(
     switch(distribution) {
     case 0  :
       for (unsigned int i = 0; i < alpha.n_slices; i++) {
-        double simsignal = alpha(0, t, i); // D and xbeta always zero
+        double simsignal = alpha(0, t, i);
         weights(i) = -0.5 * (simsignal + std::pow(y(t) / phi, 2.0) * std::exp(-simsignal)) +
           0.5 * std::pow((approx_model.y(t) - simsignal) / approx_model.H(t), 2.0);
       }
@@ -787,6 +787,9 @@ arma::mat ssm_ung::sample_model(const unsigned int predict_type) {
 
       switch(distribution) {
       case 0:
+        for (unsigned int t = 0; t < n; t++) {
+          y(0, t) = phi * exp(0.5 * alpha(0, t)) * normal(engine);
+        }
         break;
       case 1:
         for (unsigned int t = 0; t < n; t++) {
@@ -837,7 +840,6 @@ arma::cube ssm_ung::predict_past(const arma::mat& theta_posterior,
 
     update_model(theta_posterior.col(i), update_fn);
     arma::mat y(1, n);
-
     switch(distribution) {
       case 0:
         y.zeros();
@@ -873,6 +875,9 @@ arma::cube ssm_ung::predict_past(const arma::mat& theta_posterior,
 
       switch(distribution) {
       case 0:
+        for (unsigned int t = 0; t < n; t++) {
+          y(0, t) = phi * exp(0.5 * alpha(0, t, i)) * normal(engine);
+        }
         break;
       case 1:
         for (unsigned int t = 0; t < n; t++) {