From 9a5b64ca79a86998fccd006ec85c83a1bd95e190 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= <wozniak.tom@pm.me>
Date: Thu, 6 Oct 2022 10:28:32 +1100
Subject: [PATCH] cpp functions done! #15

need to be checked
---
 src/forecast.cpp | 154 +++++++++++++++++++++++++++++++++++++++++++++++
 src/forecast.h   |  37 ++++++++++++
 2 files changed, 191 insertions(+)
 create mode 100644 src/forecast.cpp
 create mode 100644 src/forecast.h

diff --git a/src/forecast.cpp b/src/forecast.cpp
new file mode 100644
index 0000000..aa22333
--- /dev/null
+++ b/src/forecast.cpp
@@ -0,0 +1,154 @@
+
+#include <RcppArmadillo.h>
+#include "sample.h"
+
+using namespace Rcpp;
+using namespace arma;
+
+
+// [[Rcpp::interfaces(cpp)]]
+// [[Rcpp::export]]
+Rcpp::List forecast_bsvar (
+    arma::cube&   posterior_B,        // (N, N, S)
+    arma::cube&   posterior_A,        // (N, K, S)
+    arma::vec&    X_T,                // (K)
+    const int     horizon
+) {
+  
+  const int       N = posterior_B.n_rows;
+  const int       S = posterior_B.n_slices;
+  const int       K = posterior_A.n_cols;
+  const int       p = (K - 1) / N;
+  
+  cube            forecasts(N, horizon, S);
+  vec             one(1, fill::value(1));
+  
+  for (int s=0; s<S; s++) {
+    
+    vec x_t       = X_T.rows(0, K - 2);
+    vec Xt        = join_cols(x_t, one);
+    
+    mat Sigma     = inv_sympd(trans(posterior_B.slice(s)) * posterior_B.slice(s));
+    
+    for (int h=1; h<horizon + 1; h++) {
+      forecasts.slice(s).col(horizon) = mvnrnd( posterior_A.slice(s) * Xt, Sigma );
+      Xt          = join_cols(forecasts.slice(s).col(horizon), Xt.rows(N, K-2), one);
+      
+    } // END h loop
+  } // END s loop
+  
+  return List::create(
+    _["forecasts"]        = forecasts
+  );
+} // END forecast_bsvar
+
+
+
+
+// [[Rcpp::interfaces(cpp)]]
+// [[Rcpp::export]]
+Rcpp::List forecast_bsvar_msh (
+    arma::cube&   posterior_B,        // (N, N, S)
+    arma::cube&   posterior_A,        // (N, K, S)
+    arma::cube&   posterior_sigma2,   // (N, M, S)
+    arma::cube&   posterior_PR_TR,    // (M, M, S)
+    arma::vec&    X_T,                // (K)
+    arma::vec&    S_T,                // (M)
+    const int     horizon
+) {
+  
+  const int       N = posterior_B.n_rows;
+  const int       M = posterior_PR_TR.n_rows;
+  const int       S = posterior_B.n_slices;
+  const int       K = posterior_A.n_cols;
+  const int       p = (K - 1) / N;
+  
+  cube            forecasts(N, horizon, S);
+  cube            forecasts_sigma2(N, horizon, S);
+  vec             one(1, fill::value(1));
+  
+  for (int s=0; s<S; s++) {
+    
+    vec x_t       = X_T.rows(0, K - 2);
+    vec Xt        = join_cols(x_t, one);
+    
+    mat B_inv     = inv(posterior_B.slice(s));
+    vec PR_ST     = S_T;
+    NumericVector zeroM   = wrap(seq_len(M) - 1);
+    
+    int St(M);
+    mat Sigma(N, N);
+    
+    for (int h=1; h<horizon + 1; h++) {
+      
+      PR_ST       = trans(posterior_PR_TR.slice(s)) * PR_ST;
+      St          = csample_num1(zeroM, wrap(PR_ST));
+      forecasts_sigma2.slice(s).col(horizon) = posterior_sigma2.slice(s).col(St);
+      
+      Sigma       = B_inv * diagmat(forecasts_sigma2.slice(s).col(horizon)) * B_inv.t();
+      forecasts.slice(s).col(horizon) = mvnrnd( posterior_A.slice(s) * Xt, Sigma );
+      Xt          = join_cols(forecasts.slice(s).col(horizon), Xt.rows(N, K-2), one);
+      
+    } // END h loop
+  } // END s loop
+  
+  return List::create(
+    _["forecasts"]        = forecasts,
+    _["forecasts_sigma"]  = sqrt(forecasts_sigma2)
+  );
+} // END forecast_bsvar_msh
+
+
+
+
+// [[Rcpp::interfaces(cpp)]]
+// [[Rcpp::export]]
+Rcpp::List forecast_bsvar_sv (
+    arma::cube&   posterior_B,        // (N, N, S)
+    arma::cube&   posterior_A,        // (N, K, S)
+    arma::cube&   posterior_h_T,      // NxS
+    arma::mat&    posterior_rho,      // NxS
+    arma::mat&    posterior_omega,    // NxS
+    arma::vec&    X_T,                // (K)
+    const int     horizon
+) {
+  
+  const int       N = posterior_B.n_rows;
+  const int       S = posterior_B.n_slices;
+  const int       K = posterior_A.n_cols;
+  const int       p = (K - 1) / N;
+  
+  cube            forecasts(N, horizon, S);
+  cube            forecasts_sigma2(N, horizon, S);
+  vec             one(1, fill::value(1));
+  
+  for (int s=0; s<S; s++) {
+    
+    vec x_t       = X_T.rows(0, K - 2);
+    vec Xt        = join_cols(x_t, one);
+    vec ht        = posterior_h_T;
+    
+    mat B_inv     = inv(posterior_B.slice(s));
+    mat Sigma(N, N);
+    double xx = 0;
+    
+    for (int h=1; h<horizon + 1; h++) {
+      
+      for (int n=0; n<N; n++) {
+        xx        = randn();
+        ht(n)     = posterior_omega(n, s) * (posterior_rho(n, s) * ht(n) + xx);
+        forecasts_sigma2(n, horizon, s) = exp(ht(n));
+      }
+      
+      Sigma       = B_inv * diagmat(forecasts_sigma2.slice(s).col(horizon)) * B_inv.t();
+      forecasts.slice(s).col(horizon) = mvnrnd( posterior_A.slice(s) * Xt, Sigma );
+      Xt          = join_cols(forecasts.slice(s).col(horizon), Xt.rows(N, K-2), one);
+      
+    } // END h loop
+  } // END s loop
+  
+  return List::create(
+    _["forecasts"]        = forecasts,
+    _["forecasts_sigma"]  = sqrt(forecasts_sigma2)
+  );
+} // END forecast_bsvar_sv
diff --git a/src/forecast.h b/src/forecast.h
new file mode 100644
index 0000000..57c56f4
--- /dev/null
+++ b/src/forecast.h
@@ -0,0 +1,37 @@
+
+#ifndef _FORECAST_H_
+#define _FORECAST_H_
+
+#include <RcppArmadillo.h>
+
+Rcpp::List forecast_bsvar (
+    arma::cube&   posterior_B,        // (N, N, S)
+    arma::cube&   posterior_A,        // (N, K, S)
+    arma::vec&    X_T,                // (K)
+    const int     horizon
+);
+
+
+Rcpp::List forecast_bsvar_msh (
+    arma::cube&   posterior_B,        // (N, N, S)
+    arma::cube&   posterior_A,        // (N, K, S)
+    arma::cube&   posterior_sigma2,   // (N, M, S)
+    arma::cube&   posterior_PR_TR,    // (M, M, S)
+    arma::vec&    X_T,                // (K)
+    arma::vec&    S_T,                // (M)
+    const int     horizon
+);
+
+
+Rcpp::List forecast_bsvar_sv (
+    arma::cube&   posterior_B,        // (N, N, S)
+    arma::cube&   posterior_A,        // (N, K, S)
+    arma::cube&   posterior_h_T,      // NxS
+    arma::mat&    posterior_rho,      // NxS
+    arma::mat&    posterior_omega,    // NxS
+    arma::vec&    X_T,                // (K)
+    const int     horizon
+);
+
+
+#endif  // _BSVARTOOLS_H_
\ No newline at end of file