From 14ff9a802679f353d122986bb0f9b8f19e424c3f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 28 Apr 2024 00:18:53 +1000 Subject: [PATCH 01/32] crated specify_prior_bvarPANEL #2 --- NAMESPACE | 1 + R/specify_bvarpanel.R | 135 +++++++++++++++++++++++++++ man/specify_prior_bvarPANEL.Rd | 164 +++++++++++++++++++++++++++++++++ 3 files changed, 300 insertions(+) create mode 100644 R/specify_bvarpanel.R create mode 100644 man/specify_prior_bvarPANEL.Rd diff --git a/NAMESPACE b/NAMESPACE index 0708706..f16da91 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,6 @@ # Generated by roxygen2: do not edit by hand +export(specify_prior_bvarPANEL) importFrom(R6,R6Class) importFrom(Rcpp,sourceCpp) useDynLib(bvarPANELs, .registration = TRUE) diff --git a/R/specify_bvarpanel.R b/R/specify_bvarpanel.R new file mode 100644 index 0000000..d7352a5 --- /dev/null +++ b/R/specify_bvarpanel.R @@ -0,0 +1,135 @@ + +#' R6 Class Representing PriorBVARPANEL +#' +#' @description +#' The class PriorBVARPANEL presents a prior specification for the Bayesian +#' hierarchical panel VAR model. +#' +#' @examples +#' prior = specify_prior_bvarPANEL$new(C = 2, N = 3, p = 1) +#' prior$M +#' +#' @export +specify_prior_bvarPANEL = R6::R6Class( + "PriorBVARPANEL", + + public = list( + + #' @field M an \code{KxN} matrix, the mean of the second-level MNIW prior + #' distribution for the global parameter matrices \eqn{\undelline{\mathbf{A}}} + #' and \eqn{\undelline{\mathbf{V}}} + M = matrix(), + + #' @field W a \code{KxK} column-specific covariance matrix of the second-level + #' MNIW prior distribution for the global parameter matrices \eqn{\undelline{\mathbf{A}}} + #' and \eqn{\undelline{\mathbf{V}}} + W = matrix(), + + #' @field S_inv an \code{NxN} row-specific precision matrix of the second-level + #' MNIW prior distribution for the global parameter matrices \eqn{\undelline{\mathbf{A}}} + #' and \eqn{\undelline{\mathbf{V}}} + S_inv = matrix(), + + #' @field S_Sigma_inv an \code{NxN} precision matrix of the second-level + #' Wishart prior distribution for the global parameter matrix \eqn{\undelline{\mathbf{\Sigma}}}. + S_Sigma_inv = matrix(), + + #' @field eta a positive shape parameter of the second-level MNIW prior distribution + #' for the global parameter matrices \eqn{\undelline{\mathbf{A}}} + #' and \eqn{\undelline{\mathbf{V}}} + eta = NA, + + #' @field mu a positive shape parameter of the second-level Wishart prior distribution + #' for the global parameter matrix \eqn{\undelline{\mathbf{\Sigma}}}. + mu = NA, + + #' @field lambda a positive shape of the second-level exp prior distribution + #' for the shape parameter \eqn{\undelline{\nu}}. + lambda = NA, + + #' @field mu_m a scalar mean of the third-level normal prior distribution + #' for the global average persistence parameter \eqn{\undelline{m}}. + mu_m = NA, + + #' @field sigma2_m a positive scalar variance of the third-level normal prior distribution + #' for the global average persistence parameter \eqn{\undelline{m}}. + sigma2_m = NA, + + #' @field s_w a positive scalar scale of the third-level gamma prior + #' distribution for parameter \eqn{\underline{w}}. + s_w = NA, + + #' @field a_w a positive scalar shape of the third-level gamma prior + #' distribution for parameter \eqn{\underline{w}}. + a_w = NA, + + #' @field s_s a positive scalar scale parameter of the third-level + #' inverted-gamma 2 prior distribution for parameter \eqn{\underline{s}}. + s_s = NA, + + #' @field nu_s a positive scalar shape parameter of the third-level + #' inverted-gamma 2 prior distribution for parameter \eqn{\underline{s}}. + nu_s = NA, + + #' @description + #' Create a new prior specification PriorBVARPANEL. + #' @param C a positive integer - the number of countries in the data. + #' @param N a positive integer - the number of dependent variables in the model. + #' @param p a positive integer - the autoregressive lag order of the SVAR model. + #' @param d a positive integer - the number of \code{exogenous} variables in the model. + #' @return A new prior specification PriorBVARPANEL. + #' @examples + #' # a prior for 2-country, 3-variable example with one lag and stationary data + #' prior = specify_prior_bvarPANEL$new(C = 2, N = 3, p = 1) + #' prior$M + #' + initialize = function(C, N, p, d = 0){ + stopifnot("Argument C must be a positive integer number." = C > 0 & C %% 1 == 0) + stopifnot("Argument N must be a positive integer number." = N > 0 & N %% 1 == 0) + stopifnot("Argument p must be a positive integer number." = p > 0 & p %% 1 == 0) + stopifnot("Argument d must be a non-negative integer number." = d >= 0 & d %% 1 == 0) + + K = N * p + 1 + d + self$M = cbind(diag(N), matrix(0, N, K - N)) + self$W = diag(c(kronecker((1:p)^2, rep(1, N) ), rep(10, d + 1))) + self$S_inv = diag(N) + self$S_Sigma_inv = diag(N) + self$eta = N + 1 + self$mu = N + 1 + self$lambda = 0.1 + self$mu_m = 1 + self$sigma2_m = 1 + self$s_w = 1 + self$a_w = 1 + self$s_s = 1 + self$nu_s = 3 + }, # END initialize + + #' @description + #' Returns the elements of the prior specification PriorBSVAR as a \code{list}. + #' + #' @examples + #' # a prior for 2-coutnry, 3-variable example with four lags + #' prior = specify_prior_bvarPANEL$new(C = 2, N = 3, p = 4) + #' prior$get_prior() # show the prior as list + #' + get_prior = function(){ + list( + M = self$M, + W = self$W, + S_inv = self$S_inv, + S_Sigma_inv = self$S_Sigma_inv, + eta = self$eta, + mu = self$mu, + lambda = self$lambda, + mu_m = self$mu_m, + sigma2_m = self$sigma2_m, + s_w = self$s_w, + a_w = self$a_w, + s_s = self$s_s, + nu_s = self$nu_s + ) + } # END get_prior + + ) # END public +) # END specify_prior_bvarPANEL diff --git a/man/specify_prior_bvarPANEL.Rd b/man/specify_prior_bvarPANEL.Rd new file mode 100644 index 0000000..b63629c --- /dev/null +++ b/man/specify_prior_bvarPANEL.Rd @@ -0,0 +1,164 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/specify_bvarpanel.R +\name{specify_prior_bvarPANEL} +\alias{specify_prior_bvarPANEL} +\title{R6 Class Representing PriorBVARPANEL} +\description{ +The class PriorBVARPANEL presents a prior specification for the Bayesian +hierarchical panel VAR model. +} +\examples{ +prior = specify_prior_bvarPANEL$new(C = 2, N = 3, p = 1) +prior$M + + +## ------------------------------------------------ +## Method `specify_prior_bvarPANEL$new` +## ------------------------------------------------ + +# a prior for 2-country, 3-variable example with one lag and stationary data +prior = specify_prior_bvarPANEL$new(C = 2, N = 3, p = 1) +prior$M + + +## ------------------------------------------------ +## Method `specify_prior_bvarPANEL$get_prior` +## ------------------------------------------------ + +# a prior for 2-coutnry, 3-variable example with four lags +prior = specify_prior_bvarPANEL$new(C = 2, N = 3, p = 4) +prior$get_prior() # show the prior as list + +} +\section{Public fields}{ +\if{html}{\out{
}} +\describe{ +\item{\code{M}}{an \code{KxN} matrix, the mean of the second-level MNIW prior +distribution for the global parameter matrices \eqn{\undelline{\mathbf{A}}} +and \eqn{\undelline{\mathbf{V}}}} + +\item{\code{W}}{a \code{KxK} column-specific covariance matrix of the second-level +MNIW prior distribution for the global parameter matrices \eqn{\undelline{\mathbf{A}}} +and \eqn{\undelline{\mathbf{V}}}} + +\item{\code{S_inv}}{an \code{NxN} row-specific precision matrix of the second-level +MNIW prior distribution for the global parameter matrices \eqn{\undelline{\mathbf{A}}} +and \eqn{\undelline{\mathbf{V}}}} + +\item{\code{S_Sigma_inv}}{an \code{NxN} precision matrix of the second-level +Wishart prior distribution for the global parameter matrix \eqn{\undelline{\mathbf{\Sigma}}}.} + +\item{\code{eta}}{a positive shape parameter of the second-level MNIW prior distribution +for the global parameter matrices \eqn{\undelline{\mathbf{A}}} +and \eqn{\undelline{\mathbf{V}}}} + +\item{\code{mu}}{a positive shape parameter of the second-level Wishart prior distribution +for the global parameter matrix \eqn{\undelline{\mathbf{\Sigma}}}.} + +\item{\code{lambda}}{a positive shape of the second-level exp prior distribution +for the shape parameter \eqn{\undelline{\nu}}.} + +\item{\code{mu_m}}{a scalar mean of the third-level normal prior distribution +for the global average persistence parameter \eqn{\undelline{m}}.} + +\item{\code{sigma2_m}}{a positive scalar variance of the third-level normal prior distribution +for the global average persistence parameter \eqn{\undelline{m}}.} + +\item{\code{s_w}}{a positive scalar scale of the third-level gamma prior +distribution for parameter \eqn{\underline{w}}.} + +\item{\code{a_w}}{a positive scalar shape of the third-level gamma prior +distribution for parameter \eqn{\underline{w}}.} + +\item{\code{s_s}}{a positive scalar scale parameter of the third-level +inverted-gamma 2 prior distribution for parameter \eqn{\underline{s}}.} + +\item{\code{nu_s}}{a positive scalar shape parameter of the third-level +inverted-gamma 2 prior distribution for parameter \eqn{\underline{s}}.} +} +\if{html}{\out{
}} +} +\section{Methods}{ +\subsection{Public methods}{ +\itemize{ +\item \href{#method-PriorBVARPANEL-new}{\code{specify_prior_bvarPANEL$new()}} +\item \href{#method-PriorBVARPANEL-get_prior}{\code{specify_prior_bvarPANEL$get_prior()}} +\item \href{#method-PriorBVARPANEL-clone}{\code{specify_prior_bvarPANEL$clone()}} +} +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-PriorBVARPANEL-new}{}}} +\subsection{Method \code{new()}}{ +Create a new prior specification PriorBVARPANEL. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_prior_bvarPANEL$new(C, N, p, d = 0)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{C}}{a positive integer - the number of countries in the data.} + +\item{\code{N}}{a positive integer - the number of dependent variables in the model.} + +\item{\code{p}}{a positive integer - the autoregressive lag order of the SVAR model.} + +\item{\code{d}}{a positive integer - the number of \code{exogenous} variables in the model.} +} +\if{html}{\out{
}} +} +\subsection{Returns}{ +A new prior specification PriorBVARPANEL. +} +\subsection{Examples}{ +\if{html}{\out{
}} +\preformatted{# a prior for 2-country, 3-variable example with one lag and stationary data +prior = specify_prior_bvarPANEL$new(C = 2, N = 3, p = 1) +prior$M + +} +\if{html}{\out{
}} + +} + +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-PriorBVARPANEL-get_prior}{}}} +\subsection{Method \code{get_prior()}}{ +Returns the elements of the prior specification PriorBSVAR as a \code{list}. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_prior_bvarPANEL$get_prior()}\if{html}{\out{
}} +} + +\subsection{Examples}{ +\if{html}{\out{
}} +\preformatted{# a prior for 2-coutnry, 3-variable example with four lags +prior = specify_prior_bvarPANEL$new(C = 2, N = 3, p = 4) +prior$get_prior() # show the prior as list + +} +\if{html}{\out{
}} + +} + +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-PriorBVARPANEL-clone}{}}} +\subsection{Method \code{clone()}}{ +The objects of this class are cloneable with this method. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_prior_bvarPANEL$clone(deep = FALSE)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{deep}}{Whether to make a deep clone.} +} +\if{html}{\out{
}} +} +} +} From d97e0564d8344fa59719bc18da5e2b472c3ba833 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 28 Apr 2024 19:22:01 +1000 Subject: [PATCH 02/32] created specify_starting_values_bvarPANEL #2 and improved notation for eqn for parameters for compatibility with tex --- NAMESPACE | 1 + R/specify_bvarpanel.R | 161 ++++++++++++++++++++--- man/specify_prior_bvarPANEL.Rd | 36 ++--- man/specify_starting_values_bvarPANEL.Rd | 146 ++++++++++++++++++++ 4 files changed, 305 insertions(+), 39 deletions(-) create mode 100644 man/specify_starting_values_bvarPANEL.Rd diff --git a/NAMESPACE b/NAMESPACE index f16da91..b497a94 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,7 @@ # Generated by roxygen2: do not edit by hand export(specify_prior_bvarPANEL) +export(specify_starting_values_bvarPANEL) importFrom(R6,R6Class) importFrom(Rcpp,sourceCpp) useDynLib(bvarPANELs, .registration = TRUE) diff --git a/R/specify_bvarpanel.R b/R/specify_bvarpanel.R index d7352a5..f584c66 100644 --- a/R/specify_bvarpanel.R +++ b/R/specify_bvarpanel.R @@ -16,59 +16,59 @@ specify_prior_bvarPANEL = R6::R6Class( public = list( #' @field M an \code{KxN} matrix, the mean of the second-level MNIW prior - #' distribution for the global parameter matrices \eqn{\undelline{\mathbf{A}}} - #' and \eqn{\undelline{\mathbf{V}}} + #' distribution for the global parameter matrices \eqn{\mathbf{A}} + #' and \eqn{\mathbf{V}} M = matrix(), #' @field W a \code{KxK} column-specific covariance matrix of the second-level - #' MNIW prior distribution for the global parameter matrices \eqn{\undelline{\mathbf{A}}} - #' and \eqn{\undelline{\mathbf{V}}} + #' MNIW prior distribution for the global parameter matrices \eqn{\mathbf{A}} + #' and \eqn{\mathbf{V}} W = matrix(), #' @field S_inv an \code{NxN} row-specific precision matrix of the second-level - #' MNIW prior distribution for the global parameter matrices \eqn{\undelline{\mathbf{A}}} - #' and \eqn{\undelline{\mathbf{V}}} + #' MNIW prior distribution for the global parameter matrices \eqn{\mathbf{A}} + #' and \eqn{\mathbf{V}} S_inv = matrix(), #' @field S_Sigma_inv an \code{NxN} precision matrix of the second-level - #' Wishart prior distribution for the global parameter matrix \eqn{\undelline{\mathbf{\Sigma}}}. + #' Wishart prior distribution for the global parameter matrix \eqn{\mathbf{\Sigma}}. S_Sigma_inv = matrix(), #' @field eta a positive shape parameter of the second-level MNIW prior distribution - #' for the global parameter matrices \eqn{\undelline{\mathbf{A}}} - #' and \eqn{\undelline{\mathbf{V}}} + #' for the global parameter matrices \eqn{\mathbf{A}} + #' and \eqn{\mathbf{V}} eta = NA, - #' @field mu a positive shape parameter of the second-level Wishart prior distribution - #' for the global parameter matrix \eqn{\undelline{\mathbf{\Sigma}}}. - mu = NA, + #' @field mu_Sigma a positive shape parameter of the second-level Wishart prior distribution + #' for the global parameter matrix \eqn{\mathbf{\Sigma}}. + mu_Sigma = NA, #' @field lambda a positive shape of the second-level exp prior distribution - #' for the shape parameter \eqn{\undelline{\nu}}. + #' for the shape parameter \eqn{\nu}. lambda = NA, #' @field mu_m a scalar mean of the third-level normal prior distribution - #' for the global average persistence parameter \eqn{\undelline{m}}. + #' for the global average persistence parameter \eqn{m}. mu_m = NA, #' @field sigma2_m a positive scalar variance of the third-level normal prior distribution - #' for the global average persistence parameter \eqn{\undelline{m}}. + #' for the global average persistence parameter \eqn{m}. sigma2_m = NA, #' @field s_w a positive scalar scale of the third-level gamma prior - #' distribution for parameter \eqn{\underline{w}}. + #' distribution for parameter \eqn{w}. s_w = NA, #' @field a_w a positive scalar shape of the third-level gamma prior - #' distribution for parameter \eqn{\underline{w}}. + #' distribution for parameter \eqn{w}. a_w = NA, #' @field s_s a positive scalar scale parameter of the third-level - #' inverted-gamma 2 prior distribution for parameter \eqn{\underline{s}}. + #' inverted-gamma 2 prior distribution for parameter \eqn{s}. s_s = NA, #' @field nu_s a positive scalar shape parameter of the third-level - #' inverted-gamma 2 prior distribution for parameter \eqn{\underline{s}}. + #' inverted-gamma 2 prior distribution for parameter \eqn{s}. nu_s = NA, #' @description @@ -95,7 +95,7 @@ specify_prior_bvarPANEL = R6::R6Class( self$S_inv = diag(N) self$S_Sigma_inv = diag(N) self$eta = N + 1 - self$mu = N + 1 + self$mu_Sigma = N + 1 self$lambda = 0.1 self$mu_m = 1 self$sigma2_m = 1 @@ -120,7 +120,7 @@ specify_prior_bvarPANEL = R6::R6Class( S_inv = self$S_inv, S_Sigma_inv = self$S_Sigma_inv, eta = self$eta, - mu = self$mu, + mu_Sigma = self$mu_Sigma, lambda = self$lambda, mu_m = self$mu_m, sigma2_m = self$sigma2_m, @@ -133,3 +133,122 @@ specify_prior_bvarPANEL = R6::R6Class( ) # END public ) # END specify_prior_bvarPANEL + + + + + + + +#' R6 Class Representing StartingValuesBVARPANEL +#' +#' @description +#' The class StartingValuesBVARPANEL presents starting values for the Bayesian +#' hierarchical panel VAR model. +#' +#' @examples +#' # starting values for a Bayesian Panel VAR +#' sv = specify_starting_values_bvarPANEL$new(C = 2, N = 3, p = 1) +#' +#' @export +specify_starting_values_bvarPANEL = R6::R6Class( + "StartingValuesBVARPANEL", + + public = list( + + #' @field A_c an \code{KxNxC} array of starting values for the local parameter + #' \eqn{\mathbf{A}_c}. + A_c = array(), + + #' @field Sigma_c an \code{NxNxC} array of starting values for the local + #' parameter \eqn{\mathbf{\Sigma}_c}. + Sigma_c = array(), + + #' @field A an \code{KxN} matrix of starting values for the global parameter + #' \eqn{\mathbf{A}}. + A = matrix(), + + #' @field V an \code{KxK} matrix of starting values for the global parameter + #' \eqn{\mathbf{V}}. + V = matrix(), + + #' @field Sigma an \code{NxN} matrix of starting values for the global parameter + #' \eqn{\mathbf{\Sigma}}. + Sigma = matrix(), + + #' @field nu a positive scalar with starting values for the global parameter + #' \eqn{\nu}. + nu = NA, + + #' @field m a positive scalar with starting values for the global hyper-parameter + #' \eqn{m}. + m = NA, + + #' @field w a positive scalar with starting values for the global hyper-parameter + #' \eqn{w}. + w = NA, + + #' @field s a positive scalar with starting values for the global hyper-parameter + #' \eqn{s}. + s = NA, + + #' @description + #' Create new starting values StartingValuesBVARPANEL + #' @param C a positive integer - the number of countries in the data. + #' @param N a positive integer - the number of dependent variables in the model. + #' @param p a positive integer - the autoregressive lag order of the SVAR model. + #' @param d a positive integer - the number of \code{exogenous} variables in the model. + #' @return Starting values StartingValuesBVARPANEL + #' @examples + #' # starting values for Bayesian Panel VAR 2-country model with 4 lags for a 3-variable system. + #' sv = specify_starting_values_bvarPANEL$new(C = 2, N = 3, p = 4) + #' + initialize = function(C, N, p, d = 0){ + stopifnot("Argument C must be a positive integer number." = C > 0 & C %% 1 == 0) + stopifnot("Argument N must be a positive integer number." = N > 0 & N %% 1 == 0) + stopifnot("Argument p must be a positive integer number." = p > 0 & p %% 1 == 0) + stopifnot("Argument d must be a non-negative integer number." = d >= 0 & d %% 1 == 0) + + K = N * p + 1 + d + self$A_c = matrix(stats::rnorm(C * K * N, sd = 0.001), c(K, N, C)) + self$Sigma_c = stats::rWishart(C, N + 1, diag(N)) + self$A = matrix(stats::rnorm(K * N, sd = 0.001), K, N) + diag(K)[,1:N] + self$V = stats::rWishart(1, K + 1, diag(K))[,,1] + self$Sigma = stats::rWishart(1, N + 1, diag(N))[,,1] + self$nu = stats::rgamma(1, 3) + self$m = stats::rnorm(1, sd = 0.001) + self$w = stats::rgamma(1, 1) + self$s = stats::rgamma(1, 1) + }, # END initialize + + #' @description + #' Returns the elements of the starting values StartingValuesBVARPANEL as a \code{list}. + #' + #' @examples + #' # starting values for a homoskedastic bsvar with 1 lag for a 3-variable system + #' sv = specify_starting_values_bvarPANEL$new(C = 2, N = 3, p = 1) + #' sv$get_starting_values() # show starting values as list + #' + get_starting_values = function(){ + list( + A_c = self$A_c, + Sigma_c = self$Sigma_c, + A = self$A, + V = self$V, + Sigma = self$Sigma, + nu = self$nu, + m = self$m, + w = self$w, + s = self$s + ) + } # END get_starting_values + ) # END public +) # END specify_starting_values_bvarPANEL + + + + + + + + diff --git a/man/specify_prior_bvarPANEL.Rd b/man/specify_prior_bvarPANEL.Rd index b63629c..75eea6b 100644 --- a/man/specify_prior_bvarPANEL.Rd +++ b/man/specify_prior_bvarPANEL.Rd @@ -34,47 +34,47 @@ prior$get_prior() # show the prior as list \if{html}{\out{
}} \describe{ \item{\code{M}}{an \code{KxN} matrix, the mean of the second-level MNIW prior -distribution for the global parameter matrices \eqn{\undelline{\mathbf{A}}} -and \eqn{\undelline{\mathbf{V}}}} +distribution for the global parameter matrices \eqn{\mathbf{A}} +and \eqn{\mathbf{V}}} \item{\code{W}}{a \code{KxK} column-specific covariance matrix of the second-level -MNIW prior distribution for the global parameter matrices \eqn{\undelline{\mathbf{A}}} -and \eqn{\undelline{\mathbf{V}}}} +MNIW prior distribution for the global parameter matrices \eqn{\mathbf{A}} +and \eqn{\mathbf{V}}} \item{\code{S_inv}}{an \code{NxN} row-specific precision matrix of the second-level -MNIW prior distribution for the global parameter matrices \eqn{\undelline{\mathbf{A}}} -and \eqn{\undelline{\mathbf{V}}}} +MNIW prior distribution for the global parameter matrices \eqn{\mathbf{A}} +and \eqn{\mathbf{V}}} \item{\code{S_Sigma_inv}}{an \code{NxN} precision matrix of the second-level -Wishart prior distribution for the global parameter matrix \eqn{\undelline{\mathbf{\Sigma}}}.} +Wishart prior distribution for the global parameter matrix \eqn{\mathbf{\Sigma}}.} \item{\code{eta}}{a positive shape parameter of the second-level MNIW prior distribution -for the global parameter matrices \eqn{\undelline{\mathbf{A}}} -and \eqn{\undelline{\mathbf{V}}}} +for the global parameter matrices \eqn{\mathbf{A}} +and \eqn{\mathbf{V}}} -\item{\code{mu}}{a positive shape parameter of the second-level Wishart prior distribution -for the global parameter matrix \eqn{\undelline{\mathbf{\Sigma}}}.} +\item{\code{mu_Sigma}}{a positive shape parameter of the second-level Wishart prior distribution +for the global parameter matrix \eqn{\mathbf{\Sigma}}.} \item{\code{lambda}}{a positive shape of the second-level exp prior distribution -for the shape parameter \eqn{\undelline{\nu}}.} +for the shape parameter \eqn{\nu}.} \item{\code{mu_m}}{a scalar mean of the third-level normal prior distribution -for the global average persistence parameter \eqn{\undelline{m}}.} +for the global average persistence parameter \eqn{m}.} \item{\code{sigma2_m}}{a positive scalar variance of the third-level normal prior distribution -for the global average persistence parameter \eqn{\undelline{m}}.} +for the global average persistence parameter \eqn{m}.} \item{\code{s_w}}{a positive scalar scale of the third-level gamma prior -distribution for parameter \eqn{\underline{w}}.} +distribution for parameter \eqn{w}.} \item{\code{a_w}}{a positive scalar shape of the third-level gamma prior -distribution for parameter \eqn{\underline{w}}.} +distribution for parameter \eqn{w}.} \item{\code{s_s}}{a positive scalar scale parameter of the third-level -inverted-gamma 2 prior distribution for parameter \eqn{\underline{s}}.} +inverted-gamma 2 prior distribution for parameter \eqn{s}.} \item{\code{nu_s}}{a positive scalar shape parameter of the third-level -inverted-gamma 2 prior distribution for parameter \eqn{\underline{s}}.} +inverted-gamma 2 prior distribution for parameter \eqn{s}.} } \if{html}{\out{
}} } diff --git a/man/specify_starting_values_bvarPANEL.Rd b/man/specify_starting_values_bvarPANEL.Rd new file mode 100644 index 0000000..b575778 --- /dev/null +++ b/man/specify_starting_values_bvarPANEL.Rd @@ -0,0 +1,146 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/specify_bvarpanel.R +\name{specify_starting_values_bvarPANEL} +\alias{specify_starting_values_bvarPANEL} +\title{R6 Class Representing StartingValuesBVARPANEL} +\description{ +The class StartingValuesBVARPANEL presents starting values for the Bayesian +hierarchical panel VAR model. +} +\examples{ +# starting values for a Bayesian Panel VAR +sv = specify_starting_values_bvarPANEL$new(C = 2, N = 3, p = 1) + + +## ------------------------------------------------ +## Method `specify_starting_values_bvarPANEL$new` +## ------------------------------------------------ + +# starting values for Bayesian Panel VAR 2-country model with 4 lags for a 3-variable system. +sv = specify_starting_values_bvarPANEL$new(C = 2, N = 3, p = 4) + + +## ------------------------------------------------ +## Method `specify_starting_values_bvarPANEL$get_starting_values` +## ------------------------------------------------ + +# starting values for a homoskedastic bsvar with 1 lag for a 3-variable system +sv = specify_starting_values_bvarPANEL$new(C = 2, N = 3, p = 1) +sv$get_starting_values() # show starting values as list + +} +\section{Public fields}{ +\if{html}{\out{
}} +\describe{ +\item{\code{A_c}}{an \code{KxNxC} array of starting values for the local parameter +\eqn{\mathbf{A}_c}.} + +\item{\code{Sigma_c}}{an \code{NxNxC} array of starting values for the local +parameter \eqn{\mathbf{\Sigma}_c}.} + +\item{\code{A}}{an \code{KxN} matrix of starting values for the global parameter +\eqn{\mathbf{A}}.} + +\item{\code{V}}{an \code{KxK} matrix of starting values for the global parameter +\eqn{\mathbf{V}}.} + +\item{\code{Sigma}}{an \code{NxN} matrix of starting values for the global parameter +\eqn{\mathbf{\Sigma}}.} + +\item{\code{nu}}{a positive scalar with starting values for the global parameter +\eqn{\nu}.} + +\item{\code{m}}{a positive scalar with starting values for the global hyper-parameter +\eqn{m}.} + +\item{\code{w}}{a positive scalar with starting values for the global hyper-parameter +\eqn{w}.} + +\item{\code{s}}{a positive scalar with starting values for the global hyper-parameter +\eqn{s}.} +} +\if{html}{\out{
}} +} +\section{Methods}{ +\subsection{Public methods}{ +\itemize{ +\item \href{#method-StartingValuesBVARPANEL-new}{\code{specify_starting_values_bvarPANEL$new()}} +\item \href{#method-StartingValuesBVARPANEL-get_starting_values}{\code{specify_starting_values_bvarPANEL$get_starting_values()}} +\item \href{#method-StartingValuesBVARPANEL-clone}{\code{specify_starting_values_bvarPANEL$clone()}} +} +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-StartingValuesBVARPANEL-new}{}}} +\subsection{Method \code{new()}}{ +Create new starting values StartingValuesBVARPANEL +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_starting_values_bvarPANEL$new(C, N, p, d = 0)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{C}}{a positive integer - the number of countries in the data.} + +\item{\code{N}}{a positive integer - the number of dependent variables in the model.} + +\item{\code{p}}{a positive integer - the autoregressive lag order of the SVAR model.} + +\item{\code{d}}{a positive integer - the number of \code{exogenous} variables in the model.} +} +\if{html}{\out{
}} +} +\subsection{Returns}{ +Starting values StartingValuesBVARPANEL +} +\subsection{Examples}{ +\if{html}{\out{
}} +\preformatted{# starting values for Bayesian Panel VAR 2-country model with 4 lags for a 3-variable system. +sv = specify_starting_values_bvarPANEL$new(C = 2, N = 3, p = 4) + +} +\if{html}{\out{
}} + +} + +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-StartingValuesBVARPANEL-get_starting_values}{}}} +\subsection{Method \code{get_starting_values()}}{ +Returns the elements of the starting values StartingValuesBVARPANEL as a \code{list}. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_starting_values_bvarPANEL$get_starting_values()}\if{html}{\out{
}} +} + +\subsection{Examples}{ +\if{html}{\out{
}} +\preformatted{# starting values for a homoskedastic bsvar with 1 lag for a 3-variable system +sv = specify_starting_values_bvarPANEL$new(C = 2, N = 3, p = 1) +sv$get_starting_values() # show starting values as list + +} +\if{html}{\out{
}} + +} + +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-StartingValuesBVARPANEL-clone}{}}} +\subsection{Method \code{clone()}}{ +The objects of this class are cloneable with this method. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_starting_values_bvarPANEL$clone(deep = FALSE)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{deep}}{Whether to make a deep clone.} +} +\if{html}{\out{
}} +} +} +} From 2dc712c791cbc05b7f73f95421109ae8aa9ac23d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 28 Apr 2024 20:18:41 +1000 Subject: [PATCH 03/32] created specify_panel_data_matrices #2 --- NAMESPACE | 1 + R/specify_bvarpanel.R | 92 +++++++++++++++++++----- man/specify_panel_data_matrices.Rd | 80 +++++++++++++++++++++ man/specify_prior_bvarPANEL.Rd | 12 ++-- man/specify_starting_values_bvarPANEL.Rd | 7 +- 5 files changed, 164 insertions(+), 28 deletions(-) create mode 100644 man/specify_panel_data_matrices.Rd diff --git a/NAMESPACE b/NAMESPACE index b497a94..e3e3d65 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,6 @@ # Generated by roxygen2: do not edit by hand +export(specify_panel_data_matrices) export(specify_prior_bvarPANEL) export(specify_starting_values_bvarPANEL) importFrom(R6,R6Class) diff --git a/R/specify_bvarpanel.R b/R/specify_bvarpanel.R index f584c66..6cce087 100644 --- a/R/specify_bvarpanel.R +++ b/R/specify_bvarpanel.R @@ -39,8 +39,8 @@ specify_prior_bvarPANEL = R6::R6Class( #' and \eqn{\mathbf{V}} eta = NA, - #' @field mu_Sigma a positive shape parameter of the second-level Wishart prior distribution - #' for the global parameter matrix \eqn{\mathbf{\Sigma}}. + #' @field mu_Sigma a positive shape parameter of the second-level Wishart prior + #' distribution for the global parameter matrix \eqn{\mathbf{\Sigma}}. mu_Sigma = NA, #' @field lambda a positive shape of the second-level exp prior distribution @@ -51,8 +51,8 @@ specify_prior_bvarPANEL = R6::R6Class( #' for the global average persistence parameter \eqn{m}. mu_m = NA, - #' @field sigma2_m a positive scalar variance of the third-level normal prior distribution - #' for the global average persistence parameter \eqn{m}. + #' @field sigma2_m a positive scalar variance of the third-level normal prior + #' distribution for the global average persistence parameter \eqn{m}. sigma2_m = NA, #' @field s_w a positive scalar scale of the third-level gamma prior @@ -76,22 +76,20 @@ specify_prior_bvarPANEL = R6::R6Class( #' @param C a positive integer - the number of countries in the data. #' @param N a positive integer - the number of dependent variables in the model. #' @param p a positive integer - the autoregressive lag order of the SVAR model. - #' @param d a positive integer - the number of \code{exogenous} variables in the model. #' @return A new prior specification PriorBVARPANEL. #' @examples #' # a prior for 2-country, 3-variable example with one lag and stationary data #' prior = specify_prior_bvarPANEL$new(C = 2, N = 3, p = 1) #' prior$M #' - initialize = function(C, N, p, d = 0){ + initialize = function(C, N, p){ stopifnot("Argument C must be a positive integer number." = C > 0 & C %% 1 == 0) stopifnot("Argument N must be a positive integer number." = N > 0 & N %% 1 == 0) stopifnot("Argument p must be a positive integer number." = p > 0 & p %% 1 == 0) - stopifnot("Argument d must be a non-negative integer number." = d >= 0 & d %% 1 == 0) - - K = N * p + 1 + d + + K = N * p + 1 self$M = cbind(diag(N), matrix(0, N, K - N)) - self$W = diag(c(kronecker((1:p)^2, rep(1, N) ), rep(10, d + 1))) + self$W = diag(c(kronecker((1:p)^2, rep(1, N) ), rep(10, 1))) self$S_inv = diag(N) self$S_Sigma_inv = diag(N) self$eta = N + 1 @@ -130,7 +128,6 @@ specify_prior_bvarPANEL = R6::R6Class( nu_s = self$nu_s ) } # END get_prior - ) # END public ) # END specify_prior_bvarPANEL @@ -197,19 +194,17 @@ specify_starting_values_bvarPANEL = R6::R6Class( #' @param C a positive integer - the number of countries in the data. #' @param N a positive integer - the number of dependent variables in the model. #' @param p a positive integer - the autoregressive lag order of the SVAR model. - #' @param d a positive integer - the number of \code{exogenous} variables in the model. #' @return Starting values StartingValuesBVARPANEL #' @examples #' # starting values for Bayesian Panel VAR 2-country model with 4 lags for a 3-variable system. #' sv = specify_starting_values_bvarPANEL$new(C = 2, N = 3, p = 4) #' - initialize = function(C, N, p, d = 0){ + initialize = function(C, N, p){ stopifnot("Argument C must be a positive integer number." = C > 0 & C %% 1 == 0) stopifnot("Argument N must be a positive integer number." = N > 0 & N %% 1 == 0) stopifnot("Argument p must be a positive integer number." = p > 0 & p %% 1 == 0) - stopifnot("Argument d must be a non-negative integer number." = d >= 0 & d %% 1 == 0) - K = N * p + 1 + d + K = N * p + 1 self$A_c = matrix(stats::rnorm(C * K * N, sd = 0.001), c(K, N, C)) self$Sigma_c = stats::rWishart(C, N + 1, diag(N)) self$A = matrix(stats::rnorm(K * N, sd = 0.001), K, N) + diag(K)[,1:N] @@ -222,7 +217,8 @@ specify_starting_values_bvarPANEL = R6::R6Class( }, # END initialize #' @description - #' Returns the elements of the starting values StartingValuesBVARPANEL as a \code{list}. + #' Returns the elements of the starting values StartingValuesBVARPANEL as + #' a \code{list}. #' #' @examples #' # starting values for a homoskedastic bsvar with 1 lag for a 3-variable system @@ -250,5 +246,67 @@ specify_starting_values_bvarPANEL = R6::R6Class( - +#' R6 Class Representing DataMatricesBVARPANEL +#' +#' @description +#' The class DataMatricesBVARPANEL presents the data matrices of dependent +#' variables, \eqn{\mathbf{Y}_c}, and regressors, \eqn{\mathbf{X}_c}, for the +#' Bayesian Panel VAR model for all countries \eqn{c = 1, ..., C}. +#' +#' @export +specify_panel_data_matrices = R6::R6Class( + "DataMatricesBVARPANEL", + + public = list( + + #' @field Y a list with \code{C} elements with \code{T_c x N} matrices of + #' dependent variables, \eqn{\mathbf{Y}_c}. + Y = list(), + + #' @field X a list with \code{C} elements with \code{T_c x K} matrices of + #' regressors, \eqn{\mathbf{X}_c}. + X = list(), + + #' @description + #' Create new data matrices DataMatricesBVARPANEL + #' @param data a list containing \code{(T_c+p)xN} matrices with country-specific + #' time series data. + #' @param p a positive integer providing model's autoregressive lag order. + #' @return New data matrices DataMatricesBVARPANEL + initialize = function(data, p = 1L) { + if (missing(data)) { + stop("Argument data has to be specified") + } else { + stopifnot("Argument data has to be a list of matrices." = is.list(data) & all(simplify2array(lapply(data, function(x){is.matrix(x) & is.numeric(x)})))) + stopifnot("Argument data has to contain matrices with the same number of columns." = length(unique(simplify2array(lapply(data, ncol)))) == 1) + stopifnot("Argument data cannot include missing values." = all(simplify2array(lapply(dd, function(x){!any(is.na(x))})))) + } + stopifnot("Argument p must be a positive integer number." = p > 0 & p %% 1 == 0) + + C = length(data) + for (c in 1:C) { + TT = nrow(data[[c]]) + T_c = TT - p + self$Y[[c]] = data[[c]][(p + 1):TT,] + + X = matrix(0, T_c, 0) + for (i in 1:p) { + X = cbind(X, data[[c]][(p + 1):TT - i,]) + } + X = cbind(X, rep(1, T_c)) + self$X = X + } # END c loop + names(self$Y) = names(self$X) = names(data) + }, # END initialize + + #' @description + #' Returns the data matrices DataMatricesBVARPANEL as a \code{list}. + get_data_matrices = function() { + list( + Y = self$Y, + X = self$X + ) + } # END get_data_matrices + ) # END public +) # END specify_panel_data_matrices diff --git a/man/specify_panel_data_matrices.Rd b/man/specify_panel_data_matrices.Rd new file mode 100644 index 0000000..284a246 --- /dev/null +++ b/man/specify_panel_data_matrices.Rd @@ -0,0 +1,80 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/specify_bvarpanel.R +\name{specify_panel_data_matrices} +\alias{specify_panel_data_matrices} +\title{R6 Class Representing DataMatricesBVARPANEL} +\description{ +The class DataMatricesBVARPANEL presents the data matrices of dependent +variables, \eqn{\mathbf{Y}_c}, and regressors, \eqn{\mathbf{X}_c}, for the +Bayesian Panel VAR model for all countries \eqn{c = 1, ..., C}. +} +\section{Public fields}{ +\if{html}{\out{
}} +\describe{ +\item{\code{Y}}{a list with \code{C} elements with \code{T_c x N} matrices of +dependent variables, \eqn{\mathbf{Y}_c}.} + +\item{\code{X}}{a list with \code{C} elements with \code{T_c x K} matrices of +regressors, \eqn{\mathbf{X}_c}.} +} +\if{html}{\out{
}} +} +\section{Methods}{ +\subsection{Public methods}{ +\itemize{ +\item \href{#method-DataMatricesBVARPANEL-new}{\code{specify_panel_data_matrices$new()}} +\item \href{#method-DataMatricesBVARPANEL-get_data_matrices}{\code{specify_panel_data_matrices$get_data_matrices()}} +\item \href{#method-DataMatricesBVARPANEL-clone}{\code{specify_panel_data_matrices$clone()}} +} +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-DataMatricesBVARPANEL-new}{}}} +\subsection{Method \code{new()}}{ +Create new data matrices DataMatricesBVARPANEL +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_panel_data_matrices$new(data, p = 1L)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{data}}{a list containing \code{(T_c+p)xN} matrices with country-specific +time series data.} + +\item{\code{p}}{a positive integer providing model's autoregressive lag order.} +} +\if{html}{\out{
}} +} +\subsection{Returns}{ +New data matrices DataMatricesBVARPANEL +} +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-DataMatricesBVARPANEL-get_data_matrices}{}}} +\subsection{Method \code{get_data_matrices()}}{ +Returns the data matrices DataMatricesBVARPANEL as a \code{list}. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_panel_data_matrices$get_data_matrices()}\if{html}{\out{
}} +} + +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-DataMatricesBVARPANEL-clone}{}}} +\subsection{Method \code{clone()}}{ +The objects of this class are cloneable with this method. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_panel_data_matrices$clone(deep = FALSE)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{deep}}{Whether to make a deep clone.} +} +\if{html}{\out{
}} +} +} +} diff --git a/man/specify_prior_bvarPANEL.Rd b/man/specify_prior_bvarPANEL.Rd index 75eea6b..cb03743 100644 --- a/man/specify_prior_bvarPANEL.Rd +++ b/man/specify_prior_bvarPANEL.Rd @@ -52,8 +52,8 @@ Wishart prior distribution for the global parameter matrix \eqn{\mathbf{\Sigma}} for the global parameter matrices \eqn{\mathbf{A}} and \eqn{\mathbf{V}}} -\item{\code{mu_Sigma}}{a positive shape parameter of the second-level Wishart prior distribution -for the global parameter matrix \eqn{\mathbf{\Sigma}}.} +\item{\code{mu_Sigma}}{a positive shape parameter of the second-level Wishart prior +distribution for the global parameter matrix \eqn{\mathbf{\Sigma}}.} \item{\code{lambda}}{a positive shape of the second-level exp prior distribution for the shape parameter \eqn{\nu}.} @@ -61,8 +61,8 @@ for the shape parameter \eqn{\nu}.} \item{\code{mu_m}}{a scalar mean of the third-level normal prior distribution for the global average persistence parameter \eqn{m}.} -\item{\code{sigma2_m}}{a positive scalar variance of the third-level normal prior distribution -for the global average persistence parameter \eqn{m}.} +\item{\code{sigma2_m}}{a positive scalar variance of the third-level normal prior +distribution for the global average persistence parameter \eqn{m}.} \item{\code{s_w}}{a positive scalar scale of the third-level gamma prior distribution for parameter \eqn{w}.} @@ -92,7 +92,7 @@ inverted-gamma 2 prior distribution for parameter \eqn{s}.} \subsection{Method \code{new()}}{ Create a new prior specification PriorBVARPANEL. \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{specify_prior_bvarPANEL$new(C, N, p, d = 0)}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{specify_prior_bvarPANEL$new(C, N, p)}\if{html}{\out{
}} } \subsection{Arguments}{ @@ -103,8 +103,6 @@ Create a new prior specification PriorBVARPANEL. \item{\code{N}}{a positive integer - the number of dependent variables in the model.} \item{\code{p}}{a positive integer - the autoregressive lag order of the SVAR model.} - -\item{\code{d}}{a positive integer - the number of \code{exogenous} variables in the model.} } \if{html}{\out{}} } diff --git a/man/specify_starting_values_bvarPANEL.Rd b/man/specify_starting_values_bvarPANEL.Rd index b575778..93f74eb 100644 --- a/man/specify_starting_values_bvarPANEL.Rd +++ b/man/specify_starting_values_bvarPANEL.Rd @@ -75,7 +75,7 @@ parameter \eqn{\mathbf{\Sigma}_c}.} \subsection{Method \code{new()}}{ Create new starting values StartingValuesBVARPANEL \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{specify_starting_values_bvarPANEL$new(C, N, p, d = 0)}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{specify_starting_values_bvarPANEL$new(C, N, p)}\if{html}{\out{
}} } \subsection{Arguments}{ @@ -86,8 +86,6 @@ Create new starting values StartingValuesBVARPANEL \item{\code{N}}{a positive integer - the number of dependent variables in the model.} \item{\code{p}}{a positive integer - the autoregressive lag order of the SVAR model.} - -\item{\code{d}}{a positive integer - the number of \code{exogenous} variables in the model.} } \if{html}{\out{}} } @@ -109,7 +107,8 @@ sv = specify_starting_values_bvarPANEL$new(C = 2, N = 3, p = 4) \if{html}{\out{}} \if{latex}{\out{\hypertarget{method-StartingValuesBVARPANEL-get_starting_values}{}}} \subsection{Method \code{get_starting_values()}}{ -Returns the elements of the starting values StartingValuesBVARPANEL as a \code{list}. +Returns the elements of the starting values StartingValuesBVARPANEL as +a \code{list}. \subsection{Usage}{ \if{html}{\out{
}}\preformatted{specify_starting_values_bvarPANEL$get_starting_values()}\if{html}{\out{
}} } From 32ecf4d06d1cf37b15f064a951dd3cf77972f718 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 28 Apr 2024 20:43:43 +1000 Subject: [PATCH 04/32] created specify_bvarPANEL #2 and set starting values --- NAMESPACE | 1 + R/specify_bvarpanel.R | 140 +++++++++++++++- man/specify_bvarPANEL.Rd | 197 +++++++++++++++++++++++ man/specify_starting_values_bvarPANEL.Rd | 49 ++++++ 4 files changed, 386 insertions(+), 1 deletion(-) create mode 100644 man/specify_bvarPANEL.Rd diff --git a/NAMESPACE b/NAMESPACE index e3e3d65..8c95ccc 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,6 @@ # Generated by roxygen2: do not edit by hand +export(specify_bvarPANEL) export(specify_panel_data_matrices) export(specify_prior_bvarPANEL) export(specify_starting_values_bvarPANEL) diff --git a/R/specify_bvarpanel.R b/R/specify_bvarpanel.R index 6cce087..36e7d54 100644 --- a/R/specify_bvarpanel.R +++ b/R/specify_bvarpanel.R @@ -237,7 +237,34 @@ specify_starting_values_bvarPANEL = R6::R6Class( w = self$w, s = self$s ) - } # END get_starting_values + }, # END get_starting_values + + #' @description + #' Returns the elements of the starting values StartingValuesBVARPANEL as a \code{list}. + #' @param last_draw a list containing the same elements as object StartingValuesBVARPANEL. + #' @return An object of class StartingValuesBVARPANEL including the last draw + #' of the current MCMC as the starting value to be passed to the continuation + #' of the MCMC estimation. + #' + #' @examples + #' sv = specify_starting_values_bvarPANEL$new(C = 2, N = 3, p = 1) + #' + #' # Modify the starting values by: + #' sv_list = sv$get_starting_values() # getting them as list + #' sv_list$A <- matrix(rnorm(12), 3, 4) # modifying the entry + #' sv$set_starting_values(sv_list) # providing to the class object + #' + set_starting_values = function(last_draw) { + self$A_c = last_draw$A_c + self$Sigma_c = last_draw$Sigma_c + self$A = last_draw$A + self$V = last_draw$V + self$Sigma = last_draw$Sigma + self$nu = last_draw$nu + self$m = last_draw$m + self$w = last_draw$w + self$s = last_draw$s + } # END set_starting_values ) # END public ) # END specify_starting_values_bvarPANEL @@ -310,3 +337,114 @@ specify_panel_data_matrices = R6::R6Class( ) # END public ) # END specify_panel_data_matrices + + + + + + +#' R6 Class representing the specification of the BVARPANEL model +#' +#' @description +#' The class BVARPANEL presents complete specification for the Bayesian Panel +#' Vector Autoregression. +#' +#' @examples +#' \dontrun{ +#' data(us_fiscal_lsuw) +#' spec = specify_bvarPANEL$new( +#' data = us_fiscal_lsuw, +#' p = 4 +#' ) +#' } +#' +#' @export +specify_bvarPANEL = R6::R6Class( + "BVARPANEL", + + public = list( + + #' @field p a non-negative integer specifying the autoregressive lag order of the model. + p = numeric(), + + #' @field prior an object PriorBSVAR with the prior specification. + prior = list(), + + #' @field data_matrices an object DataMatricesBVARPANEL with the data matrices. + data_matrices = list(), + + #' @field starting_values an object StartingValuesBVARPANEL with the starting values. + starting_values = list(), + + #' @description + #' Create a new specification of the Bayesian Panel VAR model BVARPANEL. + #' @param data a list with \code{C} elements of \code{(T_c+p)xN} matrices + #' with time series data. + #' @param p a positive integer providing model's autoregressive lag order. + #' @return A new complete specification for the Bayesian Panel VAR model BVARPANEL. + initialize = function( + data, + p = 1L + ) { + stopifnot("Argument data has to contain matrices with the same number of columns." = length(unique(simplify2array(lapply(data, ncol)))) == 1) + stopifnot("Argument p has to be a positive integer." = ((p %% 1) == 0 & p > 0)) + + self$p = p + C = length(data) + N = unique(simplify2array(lapply(data, ncol))) + + self$data_matrices = specify_panel_data_matrices$new(data, self$p) + self$prior = specify_prior_bvarPANEL$new(C, N, self$p) + self$starting_values = specify_starting_values_bvarPANEL$new(C, N, self$p) + }, # END initialize + + #' @description + #' Returns the data matrices as the DataMatricesBVARPANEL object. + #' + #' @examples + #' \dontrun{ + #' data(us_fiscal_lsuw) + #' spec = specify_bsvar$new( + #' data = us_fiscal_lsuw, + #' p = 4 + #' ) + #' spec$get_data_matrices() + #' } + get_data_matrices = function() { + self$data_matrices$clone() + }, # END get_data_matrices + + #' @description + #' Returns the prior specification as the PriorBVARPANEL object. + #' + #' @examples + #' \dontrun{ + #' data(us_fiscal_lsuw) + #' spec = specify_bsvar$new( + #' data = us_fiscal_lsuw, + #' p = 4 + #' ) + #' spec$get_prior() + #' } + get_prior = function() { + self$prior$clone() + }, # END get_prior + + #' @description + #' Returns the starting values as the StartingValuesBVARPANEL object. + #' + #' @examples + #' \dontrun{ + #' data(us_fiscal_lsuw) + #' spec = specify_bsvar$new( + #' data = us_fiscal_lsuw, + #' p = 4 + #' ) + #' spec$get_starting_values() + #' } + get_starting_values = function() { + self$starting_values$clone() + } # END get_starting_values + ) # END public +) # END specify_bvarPANEL + diff --git a/man/specify_bvarPANEL.Rd b/man/specify_bvarPANEL.Rd new file mode 100644 index 0000000..8d7e888 --- /dev/null +++ b/man/specify_bvarPANEL.Rd @@ -0,0 +1,197 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/specify_bvarpanel.R +\name{specify_bvarPANEL} +\alias{specify_bvarPANEL} +\title{R6 Class representing the specification of the BVARPANEL model} +\description{ +The class BVARPANEL presents complete specification for the Bayesian Panel +Vector Autoregression. +} +\examples{ +\dontrun{ +data(us_fiscal_lsuw) +spec = specify_bvarPANEL$new( + data = us_fiscal_lsuw, + p = 4 +) +} + + +## ------------------------------------------------ +## Method `specify_bvarPANEL$get_data_matrices` +## ------------------------------------------------ + +\dontrun{ +data(us_fiscal_lsuw) +spec = specify_bsvar$new( + data = us_fiscal_lsuw, + p = 4 +) +spec$get_data_matrices() +} + +## ------------------------------------------------ +## Method `specify_bvarPANEL$get_prior` +## ------------------------------------------------ + +\dontrun{ +data(us_fiscal_lsuw) +spec = specify_bsvar$new( + data = us_fiscal_lsuw, + p = 4 +) +spec$get_prior() +} + +## ------------------------------------------------ +## Method `specify_bvarPANEL$get_starting_values` +## ------------------------------------------------ + +\dontrun{ +data(us_fiscal_lsuw) +spec = specify_bsvar$new( + data = us_fiscal_lsuw, + p = 4 +) +spec$get_starting_values() +} +} +\section{Public fields}{ +\if{html}{\out{
}} +\describe{ +\item{\code{p}}{a non-negative integer specifying the autoregressive lag order of the model.} + +\item{\code{prior}}{an object PriorBSVAR with the prior specification.} + +\item{\code{data_matrices}}{an object DataMatricesBVARPANEL with the data matrices.} + +\item{\code{starting_values}}{an object StartingValuesBVARPANEL with the starting values.} +} +\if{html}{\out{
}} +} +\section{Methods}{ +\subsection{Public methods}{ +\itemize{ +\item \href{#method-BVARPANEL-new}{\code{specify_bvarPANEL$new()}} +\item \href{#method-BVARPANEL-get_data_matrices}{\code{specify_bvarPANEL$get_data_matrices()}} +\item \href{#method-BVARPANEL-get_prior}{\code{specify_bvarPANEL$get_prior()}} +\item \href{#method-BVARPANEL-get_starting_values}{\code{specify_bvarPANEL$get_starting_values()}} +\item \href{#method-BVARPANEL-clone}{\code{specify_bvarPANEL$clone()}} +} +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-BVARPANEL-new}{}}} +\subsection{Method \code{new()}}{ +Create a new specification of the Bayesian Panel VAR model BVARPANEL. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_bvarPANEL$new(data, p = 1L)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{data}}{a list with \code{C} elements of \code{(T_c+p)xN} matrices +with time series data.} + +\item{\code{p}}{a positive integer providing model's autoregressive lag order.} +} +\if{html}{\out{
}} +} +\subsection{Returns}{ +A new complete specification for the Bayesian Panel VAR model BVARPANEL. +} +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-BVARPANEL-get_data_matrices}{}}} +\subsection{Method \code{get_data_matrices()}}{ +Returns the data matrices as the DataMatricesBVARPANEL object. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_bvarPANEL$get_data_matrices()}\if{html}{\out{
}} +} + +\subsection{Examples}{ +\if{html}{\out{
}} +\preformatted{\dontrun{ +data(us_fiscal_lsuw) +spec = specify_bsvar$new( + data = us_fiscal_lsuw, + p = 4 +) +spec$get_data_matrices() +} +} +\if{html}{\out{
}} + +} + +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-BVARPANEL-get_prior}{}}} +\subsection{Method \code{get_prior()}}{ +Returns the prior specification as the PriorBVARPANEL object. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_bvarPANEL$get_prior()}\if{html}{\out{
}} +} + +\subsection{Examples}{ +\if{html}{\out{
}} +\preformatted{\dontrun{ +data(us_fiscal_lsuw) +spec = specify_bsvar$new( + data = us_fiscal_lsuw, + p = 4 +) +spec$get_prior() +} +} +\if{html}{\out{
}} + +} + +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-BVARPANEL-get_starting_values}{}}} +\subsection{Method \code{get_starting_values()}}{ +Returns the starting values as the StartingValuesBVARPANEL object. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_bvarPANEL$get_starting_values()}\if{html}{\out{
}} +} + +\subsection{Examples}{ +\if{html}{\out{
}} +\preformatted{\dontrun{ +data(us_fiscal_lsuw) +spec = specify_bsvar$new( + data = us_fiscal_lsuw, + p = 4 +) +spec$get_starting_values() +} +} +\if{html}{\out{
}} + +} + +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-BVARPANEL-clone}{}}} +\subsection{Method \code{clone()}}{ +The objects of this class are cloneable with this method. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_bvarPANEL$clone(deep = FALSE)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{deep}}{Whether to make a deep clone.} +} +\if{html}{\out{
}} +} +} +} diff --git a/man/specify_starting_values_bvarPANEL.Rd b/man/specify_starting_values_bvarPANEL.Rd index 93f74eb..8c821d9 100644 --- a/man/specify_starting_values_bvarPANEL.Rd +++ b/man/specify_starting_values_bvarPANEL.Rd @@ -28,6 +28,18 @@ sv = specify_starting_values_bvarPANEL$new(C = 2, N = 3, p = 4) sv = specify_starting_values_bvarPANEL$new(C = 2, N = 3, p = 1) sv$get_starting_values() # show starting values as list + +## ------------------------------------------------ +## Method `specify_starting_values_bvarPANEL$set_starting_values` +## ------------------------------------------------ + +sv = specify_starting_values_bvarPANEL$new(C = 2, N = 3, p = 1) + +# Modify the starting values by: +sv_list = sv$get_starting_values() # getting them as list +sv_list$A <- matrix(rnorm(12), 3, 4) # modifying the entry +sv$set_starting_values(sv_list) # providing to the class object + } \section{Public fields}{ \if{html}{\out{
}} @@ -66,6 +78,7 @@ parameter \eqn{\mathbf{\Sigma}_c}.} \itemize{ \item \href{#method-StartingValuesBVARPANEL-new}{\code{specify_starting_values_bvarPANEL$new()}} \item \href{#method-StartingValuesBVARPANEL-get_starting_values}{\code{specify_starting_values_bvarPANEL$get_starting_values()}} +\item \href{#method-StartingValuesBVARPANEL-set_starting_values}{\code{specify_starting_values_bvarPANEL$set_starting_values()}} \item \href{#method-StartingValuesBVARPANEL-clone}{\code{specify_starting_values_bvarPANEL$clone()}} } } @@ -124,6 +137,42 @@ sv$get_starting_values() # show starting values as list } +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-StartingValuesBVARPANEL-set_starting_values}{}}} +\subsection{Method \code{set_starting_values()}}{ +Returns the elements of the starting values StartingValuesBVARPANEL as a \code{list}. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_starting_values_bvarPANEL$set_starting_values(last_draw)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{last_draw}}{a list containing the same elements as object StartingValuesBVARPANEL.} +} +\if{html}{\out{
}} +} +\subsection{Returns}{ +An object of class StartingValuesBVARPANEL including the last draw +of the current MCMC as the starting value to be passed to the continuation +of the MCMC estimation. +} +\subsection{Examples}{ +\if{html}{\out{
}} +\preformatted{sv = specify_starting_values_bvarPANEL$new(C = 2, N = 3, p = 1) + +# Modify the starting values by: +sv_list = sv$get_starting_values() # getting them as list +sv_list$A <- matrix(rnorm(12), 3, 4) # modifying the entry +sv$set_starting_values(sv_list) # providing to the class object + +} +\if{html}{\out{
}} + +} + } \if{html}{\out{
}} \if{html}{\out{}} From 930c42ff0863319caa11256fae4b964efaf7204d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 28 Apr 2024 20:59:24 +1000 Subject: [PATCH 05/32] created specify_posterior_bvarPANEL #2 --- NAMESPACE | 1 + R/specify_bvarpanel.R | 97 ++++++++++++++++ man/specify_posterior_bvarPANEL.Rd | 172 +++++++++++++++++++++++++++++ 3 files changed, 270 insertions(+) create mode 100644 man/specify_posterior_bvarPANEL.Rd diff --git a/NAMESPACE b/NAMESPACE index 8c95ccc..5ed4047 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,6 +2,7 @@ export(specify_bvarPANEL) export(specify_panel_data_matrices) +export(specify_posterior_bvarPANEL) export(specify_prior_bvarPANEL) export(specify_starting_values_bvarPANEL) importFrom(R6,R6Class) diff --git a/R/specify_bvarpanel.R b/R/specify_bvarpanel.R index 36e7d54..f23e1a8 100644 --- a/R/specify_bvarpanel.R +++ b/R/specify_bvarpanel.R @@ -448,3 +448,100 @@ specify_bvarPANEL = R6::R6Class( ) # END public ) # END specify_bvarPANEL + + + +#' R6 Class Representing PosteriorBVARPANEL +#' +#' @description +#' The class PosteriorBVARPANEL contains posterior output and the specification +#' including the last MCMC draw for the Bayesian Panel VAR model. +#' Note that due to the thinning of the MCMC output the starting value in element +#' \code{last_draw} might not be equal to the last draw provided in +#' element \code{posterior}. +#' +#' @seealso \code{\link{specify_bvarPANEL}} +#' +#' @examples +#' \dontrun{ +#' # This is a function that is used within estimate() +#' data(us_fiscal_lsuw) +#' specification = specify_bsvar$new(us_fiscal_lsuw, p = 4) +#' set.seed(123) +#' estimate = estimate(specification, 50) +#' class(estimate) +#' } +#' @export +specify_posterior_bvarPANEL = R6::R6Class( + "PosteriorBVARPANEL", + + private = list( + normalised = FALSE + ), # END private + + public = list( + + #' @field last_draw an object of class BVARPANEL with the last draw of the + #' current MCMC run as the starting value to be passed to the continuation + #' of the MCMC estimation using \code{estimate()}. + last_draw = list(), + + #' @field posterior a list containing Bayesian estimation output. + posterior = list(), + + #' @description + #' Create a new posterior output PosteriorBVARPANEL. + #' @param specification_bvarPANEL an object of class BVARPANEL with the last + #' draw of the current MCMC run as the starting value. + #' @param posterior_bvarPANEL a list containing Bayesian estimation output. + #' @return A posterior output PosteriorBVARPANEL. + initialize = function(specification_bvarPANEL, posterior_bvarPANEL) { + + stopifnot("Argument specification_bvarPANEL must be of class BVARPANEL." = any(class(specification_bvarPANEL) == "BVARPANEL")) + stopifnot("Argument posterior_bvarPANEL must must contain MCMC output." = is.list(posterior_bvarPANEL) & is.array(posterior_bvarPANEL$A) & is.array(posterior_bvarPANEL$Sigma) & is.vector(posterior_bvarPANEL$w)) + + self$last_draw = specification_bvarPANEL + self$posterior = posterior_bvarPANEL + }, # END initialize + + #' @description + #' Returns a list containing Bayesian estimation output. + #' + #' @examples + #' \dontrun{ + #' data(us_fiscal_lsuw) + #' specification = specify_bsvar$new(us_fiscal_lsuw) + #' set.seed(123) + #' estimate = estimate(specification, 50) + #' estimate$get_posterior() + #' } + get_posterior = function(){ + self$posterior + }, # END get_posterior + + #' @description + #' Returns an object of class BVARPANEL with the last draw of the current + #' MCMC run as the starting value to be passed to the continuation of the + #' MCMC estimation using \code{estimate()}. + #' + #' @examples + #' \dontrun{ + #' data(us_fiscal_lsuw) + #' + #' # specify the model and set seed + #' specification = specify_bsvar$new(us_fiscal_lsuw, p = 4) + #' set.seed(123) + #' + #' # run the burn-in + #' burn_in = estimate(specification, 10) + #' + #' # estimate the model + #' posterior = estimate(burn_in, 10) + #' } + get_last_draw = function(){ + self$last_draw$clone() + } # END get_last_draw + + ) # END public +) # END specify_posterior_bvarPANEL + diff --git a/man/specify_posterior_bvarPANEL.Rd b/man/specify_posterior_bvarPANEL.Rd new file mode 100644 index 0000000..4a2ebe2 --- /dev/null +++ b/man/specify_posterior_bvarPANEL.Rd @@ -0,0 +1,172 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/specify_bvarpanel.R +\name{specify_posterior_bvarPANEL} +\alias{specify_posterior_bvarPANEL} +\title{R6 Class Representing PosteriorBVARPANEL} +\description{ +The class PosteriorBVARPANEL contains posterior output and the specification +including the last MCMC draw for the Bayesian Panel VAR model. +Note that due to the thinning of the MCMC output the starting value in element +\code{last_draw} might not be equal to the last draw provided in +element \code{posterior}. +} +\examples{ +\dontrun{ +# This is a function that is used within estimate() +data(us_fiscal_lsuw) +specification = specify_bsvar$new(us_fiscal_lsuw, p = 4) +set.seed(123) +estimate = estimate(specification, 50) +class(estimate) +} + +## ------------------------------------------------ +## Method `specify_posterior_bvarPANEL$get_posterior` +## ------------------------------------------------ + +\dontrun{ +data(us_fiscal_lsuw) +specification = specify_bsvar$new(us_fiscal_lsuw) +set.seed(123) +estimate = estimate(specification, 50) +estimate$get_posterior() +} + +## ------------------------------------------------ +## Method `specify_posterior_bvarPANEL$get_last_draw` +## ------------------------------------------------ + +\dontrun{ +data(us_fiscal_lsuw) + +# specify the model and set seed +specification = specify_bsvar$new(us_fiscal_lsuw, p = 4) +set.seed(123) + +# run the burn-in +burn_in = estimate(specification, 10) + +# estimate the model +posterior = estimate(burn_in, 10) +} +} +\seealso{ +\code{\link{specify_bvarPANEL}} +} +\section{Public fields}{ +\if{html}{\out{
}} +\describe{ +\item{\code{last_draw}}{an object of class BVARPANEL with the last draw of the +current MCMC run as the starting value to be passed to the continuation +of the MCMC estimation using \code{estimate()}.} + +\item{\code{posterior}}{a list containing Bayesian estimation output.} +} +\if{html}{\out{
}} +} +\section{Methods}{ +\subsection{Public methods}{ +\itemize{ +\item \href{#method-PosteriorBVARPANEL-new}{\code{specify_posterior_bvarPANEL$new()}} +\item \href{#method-PosteriorBVARPANEL-get_posterior}{\code{specify_posterior_bvarPANEL$get_posterior()}} +\item \href{#method-PosteriorBVARPANEL-get_last_draw}{\code{specify_posterior_bvarPANEL$get_last_draw()}} +\item \href{#method-PosteriorBVARPANEL-clone}{\code{specify_posterior_bvarPANEL$clone()}} +} +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-PosteriorBVARPANEL-new}{}}} +\subsection{Method \code{new()}}{ +Create a new posterior output PosteriorBVARPANEL. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_posterior_bvarPANEL$new(specification_bvarPANEL, posterior_bvarPANEL)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{specification_bvarPANEL}}{an object of class BVARPANEL with the last +draw of the current MCMC run as the starting value.} + +\item{\code{posterior_bvarPANEL}}{a list containing Bayesian estimation output.} +} +\if{html}{\out{
}} +} +\subsection{Returns}{ +A posterior output PosteriorBVARPANEL. +} +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-PosteriorBVARPANEL-get_posterior}{}}} +\subsection{Method \code{get_posterior()}}{ +Returns a list containing Bayesian estimation output. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_posterior_bvarPANEL$get_posterior()}\if{html}{\out{
}} +} + +\subsection{Examples}{ +\if{html}{\out{
}} +\preformatted{\dontrun{ +data(us_fiscal_lsuw) +specification = specify_bsvar$new(us_fiscal_lsuw) +set.seed(123) +estimate = estimate(specification, 50) +estimate$get_posterior() +} +} +\if{html}{\out{
}} + +} + +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-PosteriorBVARPANEL-get_last_draw}{}}} +\subsection{Method \code{get_last_draw()}}{ +Returns an object of class BVARPANEL with the last draw of the current +MCMC run as the starting value to be passed to the continuation of the +MCMC estimation using \code{estimate()}. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_posterior_bvarPANEL$get_last_draw()}\if{html}{\out{
}} +} + +\subsection{Examples}{ +\if{html}{\out{
}} +\preformatted{\dontrun{ +data(us_fiscal_lsuw) + +# specify the model and set seed +specification = specify_bsvar$new(us_fiscal_lsuw, p = 4) +set.seed(123) + +# run the burn-in +burn_in = estimate(specification, 10) + +# estimate the model +posterior = estimate(burn_in, 10) +} +} +\if{html}{\out{
}} + +} + +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-PosteriorBVARPANEL-clone}{}}} +\subsection{Method \code{clone()}}{ +The objects of this class are cloneable with this method. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_posterior_bvarPANEL$clone(deep = FALSE)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{deep}}{Whether to make a deep clone.} +} +\if{html}{\out{
}} +} +} +} From 637dc3cc25225d60a24e97c1e1be5bee1d040ece Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 28 Apr 2024 21:06:22 +1000 Subject: [PATCH 06/32] Update specify_bvarpanel.R #2 --- R/specify_bvarpanel.R | 1 + 1 file changed, 1 insertion(+) diff --git a/R/specify_bvarpanel.R b/R/specify_bvarpanel.R index f23e1a8..6a7b4a6 100644 --- a/R/specify_bvarpanel.R +++ b/R/specify_bvarpanel.R @@ -545,3 +545,4 @@ specify_posterior_bvarPANEL = R6::R6Class( ) # END public ) # END specify_posterior_bvarPANEL + From 258538548ceef80c51ce29e86a6f30e1810ae921 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 28 Apr 2024 21:12:45 +1000 Subject: [PATCH 07/32] essential correction in the rmniw1 sampler #2 --- src/sample_mniw.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sample_mniw.cpp b/src/sample_mniw.cpp index 47f214f..e677389 100644 --- a/src/sample_mniw.cpp +++ b/src/sample_mniw.cpp @@ -15,7 +15,7 @@ arma::field rmniw1( ) { mat Sigma = iwishrnd(S, nu); mat X_tmp = mat(size(A), fill::randn); - mat X = A + chol(S).t() * X_tmp * chol(V); + mat X = A + chol(V).t() * X_tmp * chol(Sigma); field out(2); out(0) = X; From b474c1430e7bc94267722bd4b1dc6255ebb0710f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 28 Apr 2024 22:37:02 +1000 Subject: [PATCH 08/32] corrected prior initialisation #2 --- R/specify_bvarpanel.R | 2 +- man/specify_bvarPANEL.Rd | 2 +- man/specify_panel_data_matrices.Rd | 2 +- man/specify_posterior_bvarPANEL.Rd | 2 +- man/specify_prior_bvarPANEL.Rd | 2 +- man/specify_starting_values_bvarPANEL.Rd | 2 +- 6 files changed, 6 insertions(+), 6 deletions(-) diff --git a/R/specify_bvarpanel.R b/R/specify_bvarpanel.R index 6a7b4a6..5c5e36d 100644 --- a/R/specify_bvarpanel.R +++ b/R/specify_bvarpanel.R @@ -88,7 +88,7 @@ specify_prior_bvarPANEL = R6::R6Class( stopifnot("Argument p must be a positive integer number." = p > 0 & p %% 1 == 0) K = N * p + 1 - self$M = cbind(diag(N), matrix(0, N, K - N)) + self$M = t(cbind(diag(N), matrix(0, N, K - N))) self$W = diag(c(kronecker((1:p)^2, rep(1, N) ), rep(10, 1))) self$S_inv = diag(N) self$S_Sigma_inv = diag(N) diff --git a/man/specify_bvarPANEL.Rd b/man/specify_bvarPANEL.Rd index 8d7e888..4d77395 100644 --- a/man/specify_bvarPANEL.Rd +++ b/man/specify_bvarPANEL.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/specify_bvarpanel.R +% Please edit documentation in R/specify_bvarPANEL.R \name{specify_bvarPANEL} \alias{specify_bvarPANEL} \title{R6 Class representing the specification of the BVARPANEL model} diff --git a/man/specify_panel_data_matrices.Rd b/man/specify_panel_data_matrices.Rd index 284a246..c3242ae 100644 --- a/man/specify_panel_data_matrices.Rd +++ b/man/specify_panel_data_matrices.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/specify_bvarpanel.R +% Please edit documentation in R/specify_bvarPANEL.R \name{specify_panel_data_matrices} \alias{specify_panel_data_matrices} \title{R6 Class Representing DataMatricesBVARPANEL} diff --git a/man/specify_posterior_bvarPANEL.Rd b/man/specify_posterior_bvarPANEL.Rd index 4a2ebe2..f361cd2 100644 --- a/man/specify_posterior_bvarPANEL.Rd +++ b/man/specify_posterior_bvarPANEL.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/specify_bvarpanel.R +% Please edit documentation in R/specify_bvarPANEL.R \name{specify_posterior_bvarPANEL} \alias{specify_posterior_bvarPANEL} \title{R6 Class Representing PosteriorBVARPANEL} diff --git a/man/specify_prior_bvarPANEL.Rd b/man/specify_prior_bvarPANEL.Rd index cb03743..ac0dc8c 100644 --- a/man/specify_prior_bvarPANEL.Rd +++ b/man/specify_prior_bvarPANEL.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/specify_bvarpanel.R +% Please edit documentation in R/specify_bvarPANEL.R \name{specify_prior_bvarPANEL} \alias{specify_prior_bvarPANEL} \title{R6 Class Representing PriorBVARPANEL} diff --git a/man/specify_starting_values_bvarPANEL.Rd b/man/specify_starting_values_bvarPANEL.Rd index 8c821d9..57a868d 100644 --- a/man/specify_starting_values_bvarPANEL.Rd +++ b/man/specify_starting_values_bvarPANEL.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/specify_bvarpanel.R +% Please edit documentation in R/specify_bvarPANEL.R \name{specify_starting_values_bvarPANEL} \alias{specify_starting_values_bvarPANEL} \title{R6 Class Representing StartingValuesBVARPANEL} From 74b18d3154899c0567066740f2e8ed2ab4e4c912 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 28 Apr 2024 22:37:23 +1000 Subject: [PATCH 09/32] created sample_m #2 --- R/RcppExports.R | 4 ++++ src/RcppExports.cpp | 16 ++++++++++++++++ src/sample_mniw.cpp | 30 ++++++++++++++++++++++++++++++ 3 files changed, 50 insertions(+) diff --git a/R/RcppExports.R b/R/RcppExports.R index cade656..5d31b37 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -5,3 +5,7 @@ rmniw1 <- function(A, V, S, nu) { .Call(`_bvarPANELs_rmniw1`, A, V, S, nu) } +sample_m <- function(aux_A, aux_V, aux_s, aux_w, aux_prior) { + .Call(`_bvarPANELs_sample_m`, aux_A, aux_V, aux_s, aux_w, aux_prior) +} + diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index a6c932a..8f03c67 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -25,9 +25,25 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// sample_m +double sample_m(const arma::mat& aux_A, const arma::mat& aux_V, const double& aux_s, const double& aux_w, const Rcpp::List& aux_prior); +RcppExport SEXP _bvarPANELs_sample_m(SEXP aux_ASEXP, SEXP aux_VSEXP, SEXP aux_sSEXP, SEXP aux_wSEXP, SEXP aux_priorSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type aux_A(aux_ASEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type aux_V(aux_VSEXP); + Rcpp::traits::input_parameter< const double& >::type aux_s(aux_sSEXP); + Rcpp::traits::input_parameter< const double& >::type aux_w(aux_wSEXP); + Rcpp::traits::input_parameter< const Rcpp::List& >::type aux_prior(aux_priorSEXP); + rcpp_result_gen = Rcpp::wrap(sample_m(aux_A, aux_V, aux_s, aux_w, aux_prior)); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { {"_bvarPANELs_rmniw1", (DL_FUNC) &_bvarPANELs_rmniw1, 4}, + {"_bvarPANELs_sample_m", (DL_FUNC) &_bvarPANELs_sample_m, 5}, {NULL, NULL, 0} }; diff --git a/src/sample_mniw.cpp b/src/sample_mniw.cpp index e677389..cb686c4 100644 --- a/src/sample_mniw.cpp +++ b/src/sample_mniw.cpp @@ -22,3 +22,33 @@ arma::field rmniw1( out(1) = Sigma; return out; } // END rmniw1 + + + +// [[Rcpp::export]] +double sample_m ( + const arma::mat& aux_A, // KxN + const arma::mat& aux_V, // KxK + const double& aux_s, // scalar + const double& aux_w, // scalar + const Rcpp::List& prior +) { + + const int N = aux_A.n_cols; + mat prior_S_inv = as(prior["S_inv"]); + mat prior_S = inv_sympd(prior_S_inv); + double prior_mu_m = as(prior["mu_m"]); + double prior_sigma2_m = as(prior["sigma2_m"]); + + double precision_tmp = (1 / prior_sigma2_m); + double mean_tmp = (prior_mu_m / prior_sigma2_m); + for (int n = 0; n < N; n++) { + precision_tmp += 1 / (aux_s * aux_w * prior_S(n, n) * aux_V(n, n)); + mean_tmp += aux_A(n, n) / (aux_s * aux_w * prior_S(n, n) * aux_V(n, n)); + } + double sigma2_m_bar = 1 / precision_tmp; + double mu_m_bar = sigma2_m_bar * mean_tmp; + + double out = randn( distr_param(mu_m_bar, pow(sigma2_m_bar, 0.5)) ); + return out; +} // END sample_m From 7be4dd17aabb7b66da02ff5700db7e0eeaba30b7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 28 Apr 2024 22:40:02 +1000 Subject: [PATCH 10/32] corrected sample_m #2 --- R/RcppExports.R | 4 ++-- src/RcppExports.cpp | 8 ++++---- src/sample_mniw.cpp | 2 +- src/sample_mniw.h | 9 +++++++++ 4 files changed, 16 insertions(+), 7 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 5d31b37..24f7c94 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -5,7 +5,7 @@ rmniw1 <- function(A, V, S, nu) { .Call(`_bvarPANELs_rmniw1`, A, V, S, nu) } -sample_m <- function(aux_A, aux_V, aux_s, aux_w, aux_prior) { - .Call(`_bvarPANELs_sample_m`, aux_A, aux_V, aux_s, aux_w, aux_prior) +sample_m <- function(aux_A, aux_V, aux_s, aux_w, prior) { + .Call(`_bvarPANELs_sample_m`, aux_A, aux_V, aux_s, aux_w, prior) } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 8f03c67..19dcbbd 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -26,8 +26,8 @@ BEGIN_RCPP END_RCPP } // sample_m -double sample_m(const arma::mat& aux_A, const arma::mat& aux_V, const double& aux_s, const double& aux_w, const Rcpp::List& aux_prior); -RcppExport SEXP _bvarPANELs_sample_m(SEXP aux_ASEXP, SEXP aux_VSEXP, SEXP aux_sSEXP, SEXP aux_wSEXP, SEXP aux_priorSEXP) { +double sample_m(const arma::mat& aux_A, const arma::mat& aux_V, const double& aux_s, const double& aux_w, const Rcpp::List& prior); +RcppExport SEXP _bvarPANELs_sample_m(SEXP aux_ASEXP, SEXP aux_VSEXP, SEXP aux_sSEXP, SEXP aux_wSEXP, SEXP priorSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -35,8 +35,8 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const arma::mat& >::type aux_V(aux_VSEXP); Rcpp::traits::input_parameter< const double& >::type aux_s(aux_sSEXP); Rcpp::traits::input_parameter< const double& >::type aux_w(aux_wSEXP); - Rcpp::traits::input_parameter< const Rcpp::List& >::type aux_prior(aux_priorSEXP); - rcpp_result_gen = Rcpp::wrap(sample_m(aux_A, aux_V, aux_s, aux_w, aux_prior)); + Rcpp::traits::input_parameter< const Rcpp::List& >::type prior(priorSEXP); + rcpp_result_gen = Rcpp::wrap(sample_m(aux_A, aux_V, aux_s, aux_w, prior)); return rcpp_result_gen; END_RCPP } diff --git a/src/sample_mniw.cpp b/src/sample_mniw.cpp index cb686c4..8421a6d 100644 --- a/src/sample_mniw.cpp +++ b/src/sample_mniw.cpp @@ -24,7 +24,7 @@ arma::field rmniw1( } // END rmniw1 - +// [[Rcpp:interface(cpp)]] // [[Rcpp::export]] double sample_m ( const arma::mat& aux_A, // KxN diff --git a/src/sample_mniw.h b/src/sample_mniw.h index d2d594a..1f6a057 100644 --- a/src/sample_mniw.h +++ b/src/sample_mniw.h @@ -15,4 +15,13 @@ arma::field rmniw1( ); +double sample_m ( + const arma::mat& aux_A, // KxN + const arma::mat& aux_V, // KxK + const double& aux_s, // scalar + const double& aux_w, // scalar + const Rcpp::List& prior +); + + #endif // _SAMPLE_MNIW_H_ From 0ae29186f906039be7c53bd354e80f699f6f8f1a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 28 Apr 2024 22:54:05 +1000 Subject: [PATCH 11/32] created sample_w #2 --- R/RcppExports.R | 4 ++++ src/RcppExports.cpp | 13 +++++++++++++ src/sample_mniw.cpp | 22 ++++++++++++++++++++++ src/sample_mniw.h | 5 +++++ 4 files changed, 44 insertions(+) diff --git a/R/RcppExports.R b/R/RcppExports.R index 24f7c94..1def164 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -9,3 +9,7 @@ sample_m <- function(aux_A, aux_V, aux_s, aux_w, prior) { .Call(`_bvarPANELs_sample_m`, aux_A, aux_V, aux_s, aux_w, prior) } +sample_w <- function(aux_V, prior) { + .Call(`_bvarPANELs_sample_w`, aux_V, prior) +} + diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 19dcbbd..62a6b42 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -40,10 +40,23 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// sample_w +double sample_w(const arma::mat& aux_V, const Rcpp::List& prior); +RcppExport SEXP _bvarPANELs_sample_w(SEXP aux_VSEXP, SEXP priorSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type aux_V(aux_VSEXP); + Rcpp::traits::input_parameter< const Rcpp::List& >::type prior(priorSEXP); + rcpp_result_gen = Rcpp::wrap(sample_w(aux_V, prior)); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { {"_bvarPANELs_rmniw1", (DL_FUNC) &_bvarPANELs_rmniw1, 4}, {"_bvarPANELs_sample_m", (DL_FUNC) &_bvarPANELs_sample_m, 5}, + {"_bvarPANELs_sample_w", (DL_FUNC) &_bvarPANELs_sample_w, 2}, {NULL, NULL, 0} }; diff --git a/src/sample_mniw.cpp b/src/sample_mniw.cpp index 8421a6d..487ff5f 100644 --- a/src/sample_mniw.cpp +++ b/src/sample_mniw.cpp @@ -52,3 +52,25 @@ double sample_m ( double out = randn( distr_param(mu_m_bar, pow(sigma2_m_bar, 0.5)) ); return out; } // END sample_m + + +// [[Rcpp:interface(cpp)]] +// [[Rcpp::export]] +double sample_w ( + const arma::mat& aux_V, // KxK + const Rcpp::List& prior +) { + + int K = aux_V.n_cols; + mat prior_W = as(prior["W"]); + double prior_s_w = as(prior["s_w"]); + double prior_a_w = as(prior["a_w"]); + double prior_eta = as(prior["eta"]); + + mat aux_V_inv = inv_sympd(aux_V); + double s_w_bar = prior_s_w + 0.5 * trace(aux_V_inv * prior_W); + double a_w_bar = prior_a_w + 0.5 * K * prior_eta; + + double out = randg( distr_param(a_w_bar, s_w_bar) ); + return out; +} // END sample_w diff --git a/src/sample_mniw.h b/src/sample_mniw.h index 1f6a057..97888ef 100644 --- a/src/sample_mniw.h +++ b/src/sample_mniw.h @@ -24,4 +24,9 @@ double sample_m ( ); +double sample_w ( + const arma::mat& aux_V, // KxK + const Rcpp::List& prior +); + #endif // _SAMPLE_MNIW_H_ From 015b92da0cf65a4b46e55d5e4872d4888ac5c164 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 28 Apr 2024 23:11:38 +1000 Subject: [PATCH 12/32] created sample_s #2 --- R/RcppExports.R | 4 ++++ src/RcppExports.cpp | 16 ++++++++++++++++ src/sample_mniw.cpp | 29 +++++++++++++++++++++++++++++ src/sample_mniw.h | 10 ++++++++++ 4 files changed, 59 insertions(+) diff --git a/R/RcppExports.R b/R/RcppExports.R index 1def164..dc7f6a4 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -13,3 +13,7 @@ sample_w <- function(aux_V, prior) { .Call(`_bvarPANELs_sample_w`, aux_V, prior) } +sample_s <- function(aux_A, aux_V, aux_Sigma, aux_m, prior) { + .Call(`_bvarPANELs_sample_s`, aux_A, aux_V, aux_Sigma, aux_m, prior) +} + diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 62a6b42..2740b32 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -52,11 +52,27 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// sample_s +double sample_s(const arma::mat& aux_A, const arma::mat& aux_V, const arma::mat& aux_Sigma, const double& aux_m, const Rcpp::List& prior); +RcppExport SEXP _bvarPANELs_sample_s(SEXP aux_ASEXP, SEXP aux_VSEXP, SEXP aux_SigmaSEXP, SEXP aux_mSEXP, SEXP priorSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type aux_A(aux_ASEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type aux_V(aux_VSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type aux_Sigma(aux_SigmaSEXP); + Rcpp::traits::input_parameter< const double& >::type aux_m(aux_mSEXP); + Rcpp::traits::input_parameter< const Rcpp::List& >::type prior(priorSEXP); + rcpp_result_gen = Rcpp::wrap(sample_s(aux_A, aux_V, aux_Sigma, aux_m, prior)); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { {"_bvarPANELs_rmniw1", (DL_FUNC) &_bvarPANELs_rmniw1, 4}, {"_bvarPANELs_sample_m", (DL_FUNC) &_bvarPANELs_sample_m, 5}, {"_bvarPANELs_sample_w", (DL_FUNC) &_bvarPANELs_sample_w, 2}, + {"_bvarPANELs_sample_s", (DL_FUNC) &_bvarPANELs_sample_s, 5}, {NULL, NULL, 0} }; diff --git a/src/sample_mniw.cpp b/src/sample_mniw.cpp index 487ff5f..d8e4350 100644 --- a/src/sample_mniw.cpp +++ b/src/sample_mniw.cpp @@ -74,3 +74,32 @@ double sample_w ( double out = randg( distr_param(a_w_bar, s_w_bar) ); return out; } // END sample_w + + +// [[Rcpp:interface(cpp)]] +// [[Rcpp::export]] +double sample_s ( + const arma::mat& aux_A, // KxN + const arma::mat& aux_V, // KxK + const arma::mat& aux_Sigma, // NxN + const double& aux_m, // scalar + const Rcpp::List& prior +) { + + int N = aux_A.n_cols; + int K = aux_A.n_rows; + mat prior_M = as(prior["M"]); + mat prior_S_inv = as(prior["S_inv"]); + mat prior_S_Sigma_inv = as(prior["S_Sigma_inv"]); + double prior_s_s = as(prior["s_s"]); + double prior_nu_s = as(prior["nu_s"]); + double prior_mu_Sigma = as(prior["mu_Sigma"]); + + mat quad_tmp1 = (aux_A - aux_m * prior_M) * prior_S_inv * trans(aux_A - aux_m * prior_M); + + double s_s_bar = prior_s_s + trace(solve(aux_V, quad_tmp1)) + trace(prior_S_Sigma_inv * aux_Sigma); + double nu_s_bar = prior_nu_s + K * N + N * prior_mu_Sigma; + + double out = s_s_bar / chi2rnd( nu_s_bar ); + return out; +} // END sample_s diff --git a/src/sample_mniw.h b/src/sample_mniw.h index 97888ef..bdfa7ce 100644 --- a/src/sample_mniw.h +++ b/src/sample_mniw.h @@ -29,4 +29,14 @@ double sample_w ( const Rcpp::List& prior ); + +double sample_s ( + const arma::mat& aux_A, // KxN + const arma::mat& aux_V, // KxK + const arma::mat& aux_Sigma, // NxN + const double& aux_m, // scalar + const Rcpp::List& prior +); + + #endif // _SAMPLE_MNIW_H_ From f5b68662e246190d390ceb89ba2e975098aa72c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 28 Apr 2024 23:33:02 +1000 Subject: [PATCH 13/32] created sample_Sigma #2 --- R/RcppExports.R | 4 ++++ src/RcppExports.cpp | 15 +++++++++++++++ src/sample_mniw.cpp | 27 +++++++++++++++++++++++++++ src/sample_mniw.h | 8 ++++++++ 4 files changed, 54 insertions(+) diff --git a/R/RcppExports.R b/R/RcppExports.R index dc7f6a4..201a161 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -17,3 +17,7 @@ sample_s <- function(aux_A, aux_V, aux_Sigma, aux_m, prior) { .Call(`_bvarPANELs_sample_s`, aux_A, aux_V, aux_Sigma, aux_m, prior) } +sample_Sigma <- function(aux_Sigma_c_inv, aux_s, aux_nu, prior) { + .Call(`_bvarPANELs_sample_Sigma`, aux_Sigma_c_inv, aux_s, aux_nu, prior) +} + diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 2740b32..aad280c 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -67,12 +67,27 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// sample_Sigma +arma::mat sample_Sigma(const arma::cube& aux_Sigma_c_inv, const double& aux_s, const double& aux_nu, const Rcpp::List& prior); +RcppExport SEXP _bvarPANELs_sample_Sigma(SEXP aux_Sigma_c_invSEXP, SEXP aux_sSEXP, SEXP aux_nuSEXP, SEXP priorSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::cube& >::type aux_Sigma_c_inv(aux_Sigma_c_invSEXP); + Rcpp::traits::input_parameter< const double& >::type aux_s(aux_sSEXP); + Rcpp::traits::input_parameter< const double& >::type aux_nu(aux_nuSEXP); + Rcpp::traits::input_parameter< const Rcpp::List& >::type prior(priorSEXP); + rcpp_result_gen = Rcpp::wrap(sample_Sigma(aux_Sigma_c_inv, aux_s, aux_nu, prior)); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { {"_bvarPANELs_rmniw1", (DL_FUNC) &_bvarPANELs_rmniw1, 4}, {"_bvarPANELs_sample_m", (DL_FUNC) &_bvarPANELs_sample_m, 5}, {"_bvarPANELs_sample_w", (DL_FUNC) &_bvarPANELs_sample_w, 2}, {"_bvarPANELs_sample_s", (DL_FUNC) &_bvarPANELs_sample_s, 5}, + {"_bvarPANELs_sample_Sigma", (DL_FUNC) &_bvarPANELs_sample_Sigma, 4}, {NULL, NULL, 0} }; diff --git a/src/sample_mniw.cpp b/src/sample_mniw.cpp index d8e4350..4c5e581 100644 --- a/src/sample_mniw.cpp +++ b/src/sample_mniw.cpp @@ -103,3 +103,30 @@ double sample_s ( double out = s_s_bar / chi2rnd( nu_s_bar ); return out; } // END sample_s + + +// [[Rcpp:interface(cpp)]] +// [[Rcpp::export]] +arma::mat sample_Sigma ( + const arma::cube& aux_Sigma_c_inv, // NxNxC + const double& aux_s, // scalar + const double& aux_nu, // scalar + const Rcpp::List& prior +) { + + int C = aux_Sigma_c_inv.n_slices; + int N = aux_Sigma_c_inv.n_rows; + mat prior_S_Sigma_inv = as(prior["S_Sigma_inv"]); + double prior_mu_Sigma = as(prior["mu_Sigma"]); + + mat sum_aux_Sigma_c_inv(N, N); + for (int c = 0; c < C; c++) { + sum_aux_Sigma_c_inv += aux_Sigma_c_inv.slice(c); + } + + mat S_Sigma_bar = (prior_S_Sigma_inv / aux_s) + sum_aux_Sigma_c_inv; + double mu_bar = C * aux_nu + prior_mu_Sigma; + + mat out = wishrnd( S_Sigma_bar, mu_bar ); + return out; +} // END sample_Sigma diff --git a/src/sample_mniw.h b/src/sample_mniw.h index bdfa7ce..cdd04c7 100644 --- a/src/sample_mniw.h +++ b/src/sample_mniw.h @@ -39,4 +39,12 @@ double sample_s ( ); +arma::mat sample_Sigma ( + const arma::cube& aux_Sigma_c_inv, // NxNxC + const double& aux_s, // scalar + const double& aux_nu, // scalar + const Rcpp::List& prior +); + + #endif // _SAMPLE_MNIW_H_ From 1de8e6b431dbcde45ee35cdeb0feb2375bd4d533 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Mon, 29 Apr 2024 00:02:17 +1000 Subject: [PATCH 14/32] created sample_AV #2 --- R/RcppExports.R | 4 ++++ src/RcppExports.cpp | 17 +++++++++++++++++ src/sample_mniw.cpp | 42 ++++++++++++++++++++++++++++++++++++++++++ src/sample_mniw.h | 10 ++++++++++ 4 files changed, 73 insertions(+) diff --git a/R/RcppExports.R b/R/RcppExports.R index 201a161..0e71b04 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -21,3 +21,7 @@ sample_Sigma <- function(aux_Sigma_c_inv, aux_s, aux_nu, prior) { .Call(`_bvarPANELs_sample_Sigma`, aux_Sigma_c_inv, aux_s, aux_nu, prior) } +sample_AV <- function(aux_A_c, aux_Sigma_c_inv, aux_s, aux_m, aux_w, prior) { + .Call(`_bvarPANELs_sample_AV`, aux_A_c, aux_Sigma_c_inv, aux_s, aux_m, aux_w, prior) +} + diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index aad280c..ae59cb8 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -81,6 +81,22 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// sample_AV +arma::field sample_AV(const arma::cube& aux_A_c, const arma::cube& aux_Sigma_c_inv, const double& aux_s, const double& aux_m, const double& aux_w, const Rcpp::List& prior); +RcppExport SEXP _bvarPANELs_sample_AV(SEXP aux_A_cSEXP, SEXP aux_Sigma_c_invSEXP, SEXP aux_sSEXP, SEXP aux_mSEXP, SEXP aux_wSEXP, SEXP priorSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::cube& >::type aux_A_c(aux_A_cSEXP); + Rcpp::traits::input_parameter< const arma::cube& >::type aux_Sigma_c_inv(aux_Sigma_c_invSEXP); + Rcpp::traits::input_parameter< const double& >::type aux_s(aux_sSEXP); + Rcpp::traits::input_parameter< const double& >::type aux_m(aux_mSEXP); + Rcpp::traits::input_parameter< const double& >::type aux_w(aux_wSEXP); + Rcpp::traits::input_parameter< const Rcpp::List& >::type prior(priorSEXP); + rcpp_result_gen = Rcpp::wrap(sample_AV(aux_A_c, aux_Sigma_c_inv, aux_s, aux_m, aux_w, prior)); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { {"_bvarPANELs_rmniw1", (DL_FUNC) &_bvarPANELs_rmniw1, 4}, @@ -88,6 +104,7 @@ static const R_CallMethodDef CallEntries[] = { {"_bvarPANELs_sample_w", (DL_FUNC) &_bvarPANELs_sample_w, 2}, {"_bvarPANELs_sample_s", (DL_FUNC) &_bvarPANELs_sample_s, 5}, {"_bvarPANELs_sample_Sigma", (DL_FUNC) &_bvarPANELs_sample_Sigma, 4}, + {"_bvarPANELs_sample_AV", (DL_FUNC) &_bvarPANELs_sample_AV, 6}, {NULL, NULL, 0} }; diff --git a/src/sample_mniw.cpp b/src/sample_mniw.cpp index 4c5e581..26027c1 100644 --- a/src/sample_mniw.cpp +++ b/src/sample_mniw.cpp @@ -130,3 +130,45 @@ arma::mat sample_Sigma ( mat out = wishrnd( S_Sigma_bar, mu_bar ); return out; } // END sample_Sigma + + +// [[Rcpp:interface(cpp)]] +// [[Rcpp::export]] +arma::field sample_AV ( + const arma::cube& aux_A_c, // KxNxC + const arma::cube& aux_Sigma_c_inv, // NxNxC + const double& aux_s, // scalar + const double& aux_m, // scalar + const double& aux_w, // scalar + const Rcpp::List& prior +) { + + int C = aux_A_c.n_slices; + int N = aux_A_c.n_cols; + int K = aux_A_c.n_rows; + + mat prior_S_inv = as(prior["S_inv"]); + mat prior_M = as(prior["M"]); + mat prior_W = as(prior["W"]); + double prior_eta = as(prior["eta"]); + + mat sum_Sc_inv(N, N); + mat sum_Sc_invAt(N, K); + mat sum_ASc_invAt(K, K); + for (int c = 0; c < C; c++) { + sum_Sc_inv += aux_Sigma_c_inv.slice(c); + mat Sc_invAt = aux_Sigma_c_inv.slice(c) * aux_A_c.slice(c).t(); + sum_Sc_invAt += Sc_invAt; + sum_ASc_invAt += aux_A_c.slice(c) * Sc_invAt; + } // END c loop + + mat S_bar_inv = (prior_S_inv / aux_s) + sum_Sc_inv; + mat S_bar = inv_sympd(S_bar_inv); + mat M_bar_trans = S_bar * ( (aux_m / aux_s) * (prior_S_inv * prior_M.t()) + sum_Sc_invAt); + mat W_bar = (aux_w * prior_W) + (pow(aux_m, 2) / aux_s) * (prior_M * prior_S_inv * prior_M.t()) + + sum_ASc_invAt - M_bar_trans.t() * S_bar_inv * M_bar_trans; + double eta_bar = C * N + prior_eta; + + arma::field aux_AV = rmniw1( M_bar_trans.t(), W_bar, S_bar, eta_bar ); + return aux_AV; +} // END sample_AV diff --git a/src/sample_mniw.h b/src/sample_mniw.h index cdd04c7..876e4a1 100644 --- a/src/sample_mniw.h +++ b/src/sample_mniw.h @@ -47,4 +47,14 @@ arma::mat sample_Sigma ( ); +arma::field sample_AV ( + const arma::cube& aux_A_c, // KxNxC + const arma::cube& aux_Sigma_c_inv, // NxNxC + const double& aux_s, // scalar + const double& aux_m, // scalar + const double& aux_w, // scalar + const Rcpp::List& prior +); + + #endif // _SAMPLE_MNIW_H_ From 30553f73763257fed38ab31eac006a29683dd4f0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 5 May 2024 16:02:57 +1000 Subject: [PATCH 15/32] created sampler for A_c and Sigma_c #2 + and corrected AV sampler ESSENTIAL --- R/RcppExports.R | 4 ++++ src/RcppExports.cpp | 17 +++++++++++++++++ src/sample_mniw.cpp | 30 ++++++++++++++++++++++++++++-- src/sample_mniw.h | 10 ++++++++++ 4 files changed, 59 insertions(+), 2 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 0e71b04..2e150f8 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -25,3 +25,7 @@ sample_AV <- function(aux_A_c, aux_Sigma_c_inv, aux_s, aux_m, aux_w, prior) { .Call(`_bvarPANELs_sample_AV`, aux_A_c, aux_Sigma_c_inv, aux_s, aux_m, aux_w, prior) } +sample_A_c_Sigma_c <- function(Y_c, X_c, aux_A, aux_V, aux_Sigma, aux_nu) { + .Call(`_bvarPANELs_sample_A_c_Sigma_c`, Y_c, X_c, aux_A, aux_V, aux_Sigma, aux_nu) +} + diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index ae59cb8..675f369 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -97,6 +97,22 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// sample_A_c_Sigma_c +arma::field sample_A_c_Sigma_c(const arma::mat& Y_c, const arma::mat& X_c, const arma::mat& aux_A, const arma::mat& aux_V, const arma::mat& aux_Sigma, const double& aux_nu); +RcppExport SEXP _bvarPANELs_sample_A_c_Sigma_c(SEXP Y_cSEXP, SEXP X_cSEXP, SEXP aux_ASEXP, SEXP aux_VSEXP, SEXP aux_SigmaSEXP, SEXP aux_nuSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type Y_c(Y_cSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type X_c(X_cSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type aux_A(aux_ASEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type aux_V(aux_VSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type aux_Sigma(aux_SigmaSEXP); + Rcpp::traits::input_parameter< const double& >::type aux_nu(aux_nuSEXP); + rcpp_result_gen = Rcpp::wrap(sample_A_c_Sigma_c(Y_c, X_c, aux_A, aux_V, aux_Sigma, aux_nu)); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { {"_bvarPANELs_rmniw1", (DL_FUNC) &_bvarPANELs_rmniw1, 4}, @@ -105,6 +121,7 @@ static const R_CallMethodDef CallEntries[] = { {"_bvarPANELs_sample_s", (DL_FUNC) &_bvarPANELs_sample_s, 5}, {"_bvarPANELs_sample_Sigma", (DL_FUNC) &_bvarPANELs_sample_Sigma, 4}, {"_bvarPANELs_sample_AV", (DL_FUNC) &_bvarPANELs_sample_AV, 6}, + {"_bvarPANELs_sample_A_c_Sigma_c", (DL_FUNC) &_bvarPANELs_sample_A_c_Sigma_c, 6}, {NULL, NULL, 0} }; diff --git a/src/sample_mniw.cpp b/src/sample_mniw.cpp index 26027c1..ccf1b6a 100644 --- a/src/sample_mniw.cpp +++ b/src/sample_mniw.cpp @@ -161,7 +161,7 @@ arma::field sample_AV ( sum_Sc_invAt += Sc_invAt; sum_ASc_invAt += aux_A_c.slice(c) * Sc_invAt; } // END c loop - + mat S_bar_inv = (prior_S_inv / aux_s) + sum_Sc_inv; mat S_bar = inv_sympd(S_bar_inv); mat M_bar_trans = S_bar * ( (aux_m / aux_s) * (prior_S_inv * prior_M.t()) + sum_Sc_invAt); @@ -169,6 +169,32 @@ arma::field sample_AV ( + sum_ASc_invAt - M_bar_trans.t() * S_bar_inv * M_bar_trans; double eta_bar = C * N + prior_eta; - arma::field aux_AV = rmniw1( M_bar_trans.t(), W_bar, S_bar, eta_bar ); + arma::field aux_AV = rmniw1( M_bar_trans, S_bar, W_bar, eta_bar ); + aux_AV(0) = trans(aux_AV(0)); return aux_AV; } // END sample_AV + + + +// [[Rcpp:interface(cpp)]] +// [[Rcpp::export]] +arma::field sample_A_c_Sigma_c ( + const arma::mat& Y_c, // T_cxN + const arma::mat& X_c, // T_cxK + const arma::mat& aux_A, // KxN + const arma::mat& aux_V, // KxK + const arma::mat& aux_Sigma, // NxN + const double& aux_nu // scalar +) { + int T_c = Y_c.n_rows; + + mat aux_V_inv = inv_sympd( aux_V ); + mat V_bar_inv = X_c.t() * X_c + aux_V_inv; + mat V_bar = inv_sympd( V_bar_inv ); + mat A_bar = V_bar * ( X_c.t() * Y_c + aux_V_inv * aux_A ); + mat Sigma_bar = aux_Sigma + Y_c.t() * Y_c + aux_A.t() * aux_V_inv * aux_A - A_bar.t() * V_bar_inv * A_bar; + double nu_bar = T_c + aux_nu; + + arma::field aux_A_c_Sigma_c = rmniw1( A_bar, V_bar, Sigma_bar, nu_bar ); + return aux_A_c_Sigma_c; +} // END sample_A_c_Sigma_c diff --git a/src/sample_mniw.h b/src/sample_mniw.h index 876e4a1..3cd96eb 100644 --- a/src/sample_mniw.h +++ b/src/sample_mniw.h @@ -57,4 +57,14 @@ arma::field sample_AV ( ); +arma::field sample_A_c_Sigma_c ( + const arma::mat& Y_c, // T_cxN + const arma::mat& X_c, // T_cxK + const arma::mat& aux_A, // KxN + const arma::mat& aux_V, // KxK + const arma::mat& aux_Sigma, // NxN + const double& aux_nu // scalar +); + + #endif // _SAMPLE_MNIW_H_ From 612ed72ffb0f51a713c6eb3c29eda945060d50c5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 5 May 2024 18:13:26 +1000 Subject: [PATCH 16/32] included log_kernel_nu #2 --- R/RcppExports.R | 4 ++++ src/RcppExports.cpp | 18 ++++++++++++++++++ src/sample_mniw.cpp | 43 ++++++++++++++++++++++++++++++++++++++++++- src/sample_mniw.h | 11 +++++++++++ 4 files changed, 75 insertions(+), 1 deletion(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 2e150f8..83ee296 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -17,6 +17,10 @@ sample_s <- function(aux_A, aux_V, aux_Sigma, aux_m, prior) { .Call(`_bvarPANELs_sample_s`, aux_A, aux_V, aux_Sigma, aux_m, prior) } +log_kernel_nu <- function(aux_nu, aux_Sigma_c, aux_Sigma, prior_lambda, C, N, K) { + .Call(`_bvarPANELs_log_kernel_nu`, aux_nu, aux_Sigma_c, aux_Sigma, prior_lambda, C, N, K) +} + sample_Sigma <- function(aux_Sigma_c_inv, aux_s, aux_nu, prior) { .Call(`_bvarPANELs_sample_Sigma`, aux_Sigma_c_inv, aux_s, aux_nu, prior) } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 675f369..c95a403 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -67,6 +67,23 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// log_kernel_nu +double log_kernel_nu(const double& aux_nu, const arma::cube& aux_Sigma_c, const arma::mat& aux_Sigma, const double& prior_lambda, const int& C, const int& N, const int& K); +RcppExport SEXP _bvarPANELs_log_kernel_nu(SEXP aux_nuSEXP, SEXP aux_Sigma_cSEXP, SEXP aux_SigmaSEXP, SEXP prior_lambdaSEXP, SEXP CSEXP, SEXP NSEXP, SEXP KSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const double& >::type aux_nu(aux_nuSEXP); + Rcpp::traits::input_parameter< const arma::cube& >::type aux_Sigma_c(aux_Sigma_cSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type aux_Sigma(aux_SigmaSEXP); + Rcpp::traits::input_parameter< const double& >::type prior_lambda(prior_lambdaSEXP); + Rcpp::traits::input_parameter< const int& >::type C(CSEXP); + Rcpp::traits::input_parameter< const int& >::type N(NSEXP); + Rcpp::traits::input_parameter< const int& >::type K(KSEXP); + rcpp_result_gen = Rcpp::wrap(log_kernel_nu(aux_nu, aux_Sigma_c, aux_Sigma, prior_lambda, C, N, K)); + return rcpp_result_gen; +END_RCPP +} // sample_Sigma arma::mat sample_Sigma(const arma::cube& aux_Sigma_c_inv, const double& aux_s, const double& aux_nu, const Rcpp::List& prior); RcppExport SEXP _bvarPANELs_sample_Sigma(SEXP aux_Sigma_c_invSEXP, SEXP aux_sSEXP, SEXP aux_nuSEXP, SEXP priorSEXP) { @@ -119,6 +136,7 @@ static const R_CallMethodDef CallEntries[] = { {"_bvarPANELs_sample_m", (DL_FUNC) &_bvarPANELs_sample_m, 5}, {"_bvarPANELs_sample_w", (DL_FUNC) &_bvarPANELs_sample_w, 2}, {"_bvarPANELs_sample_s", (DL_FUNC) &_bvarPANELs_sample_s, 5}, + {"_bvarPANELs_log_kernel_nu", (DL_FUNC) &_bvarPANELs_log_kernel_nu, 7}, {"_bvarPANELs_sample_Sigma", (DL_FUNC) &_bvarPANELs_sample_Sigma, 4}, {"_bvarPANELs_sample_AV", (DL_FUNC) &_bvarPANELs_sample_AV, 6}, {"_bvarPANELs_sample_A_c_Sigma_c", (DL_FUNC) &_bvarPANELs_sample_A_c_Sigma_c, 6}, diff --git a/src/sample_mniw.cpp b/src/sample_mniw.cpp index ccf1b6a..7dc0be7 100644 --- a/src/sample_mniw.cpp +++ b/src/sample_mniw.cpp @@ -105,6 +105,48 @@ double sample_s ( } // END sample_s + +// [[Rcpp:interface(cpp)]] +// [[Rcpp::export]] +double log_kernel_nu ( + const double& aux_nu, // scalar + const arma::cube& aux_Sigma_c, // NxNxC + const arma::mat& aux_Sigma, // NxN + const double& prior_lambda, // scalar + const int& C, // scalar + const int& N, // scalar + const int& K // scalar +) { + + double log_kernel_nu = 0; + + log_kernel_nu -= 0.5 * C * N * (K + aux_nu); + log_kernel_nu -= prior_lambda * aux_nu; + double ldS = log_det_sympd(aux_Sigma); + log_kernel_nu += 0.5 * C * aux_nu * ldS; + + for (int n = 1; n < N + 1; n++) { + log_kernel_nu -= C * R::lgammafn(0.5 * (aux_nu + 1 - n)); + } // EDN n loop + + for (int c = 0; c < C; c++) { + double ldSc = log_det_sympd(aux_Sigma_c.slice(c)); + log_kernel_nu -= 0.5 * (aux_nu + N + K + 1) * ldSc; + } // END c loop + + return log_kernel_nu; +} // END log_kernel_nu + + + + + + + + + + + // [[Rcpp:interface(cpp)]] // [[Rcpp::export]] arma::mat sample_Sigma ( @@ -175,7 +217,6 @@ arma::field sample_AV ( } // END sample_AV - // [[Rcpp:interface(cpp)]] // [[Rcpp::export]] arma::field sample_A_c_Sigma_c ( diff --git a/src/sample_mniw.h b/src/sample_mniw.h index 3cd96eb..bcb437c 100644 --- a/src/sample_mniw.h +++ b/src/sample_mniw.h @@ -39,6 +39,17 @@ double sample_s ( ); +double log_kernel_nu ( + const double& aux_nu, // scalar + const arma::cube& aux_Sigma_c, // NxNxC + const arma::mat& aux_Sigma, // NxN + const double& prior_lambda, // scalar + const int& C, // scalar + const int& N, // scalar + const int& K // scalar +); + + arma::mat sample_Sigma ( const arma::cube& aux_Sigma_c_inv, // NxNxC const double& aux_s, // scalar From e10cd5c310e0e4b2cc75c54ccff261b1ba04a2f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 5 May 2024 18:25:06 +1000 Subject: [PATCH 17/32] included dependencies on RcppTN #2 to later sample nu --- DESCRIPTION | 7 ++++--- NAMESPACE | 2 ++ R/bvarPANELs-package.R | 1 + src/sample_mniw.cpp | 4 +++- 4 files changed, 10 insertions(+), 4 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 7f3d960..9d4e43e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -8,8 +8,9 @@ Author: c(person(given="Tomasz", family="Woźniak", email="wozniak.tom@pm.me", r Maintainer: Tomasz Woźniak License: GPL (>= 3) Imports: - R6, - Rcpp (>= 1.0.12) -LinkingTo: Rcpp, RcppArmadillo + Rcpp (>= 1.0.12), + RcppTN, + R6 +LinkingTo: Rcpp, RcppArmadillo, RcppTN Encoding: UTF-8 RoxygenNote: 7.3.1 diff --git a/NAMESPACE b/NAMESPACE index 5ed4047..0d3b7c0 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -7,4 +7,6 @@ export(specify_prior_bvarPANEL) export(specify_starting_values_bvarPANEL) importFrom(R6,R6Class) importFrom(Rcpp,sourceCpp) +importFrom(RcppTN,dtn) +importFrom(RcppTN,rtn) useDynLib(bvarPANELs, .registration = TRUE) diff --git a/R/bvarPANELs-package.R b/R/bvarPANELs-package.R index f918c81..936837a 100644 --- a/R/bvarPANELs-package.R +++ b/R/bvarPANELs-package.R @@ -32,6 +32,7 @@ #' @useDynLib bvarPANELs, .registration = TRUE #' @importFrom Rcpp sourceCpp #' @importFrom R6 R6Class +#' @importFrom RcppTN rtn dtn #' @note This package is currently in active development. #' @author Tomasz Woźniak \email{wozniak.tom@pm.me} #' @keywords package models ts diff --git a/src/sample_mniw.cpp b/src/sample_mniw.cpp index 7dc0be7..18f81ac 100644 --- a/src/sample_mniw.cpp +++ b/src/sample_mniw.cpp @@ -1,5 +1,5 @@ #include - +#include using namespace Rcpp; using namespace arma; @@ -147,6 +147,8 @@ double log_kernel_nu ( + + // [[Rcpp:interface(cpp)]] // [[Rcpp::export]] arma::mat sample_Sigma ( From 00e4827ad9c107707aeb9506ff379be1fa6a0c00 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 5 May 2024 18:51:25 +1000 Subject: [PATCH 18/32] created sample_nu #2 --- R/RcppExports.R | 4 ++++ src/RcppExports.cpp | 15 +++++++++++++++ src/sample_mniw.cpp | 44 +++++++++++++++++++++++++++++++++++++++----- src/sample_mniw.h | 8 ++++++++ 4 files changed, 66 insertions(+), 5 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 83ee296..8e6f27d 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -21,6 +21,10 @@ log_kernel_nu <- function(aux_nu, aux_Sigma_c, aux_Sigma, prior_lambda, C, N, K) .Call(`_bvarPANELs_log_kernel_nu`, aux_nu, aux_Sigma_c, aux_Sigma, prior_lambda, C, N, K) } +sample_nu <- function(aux_nu, aux_Sigma_c, aux_Sigma, prior) { + .Call(`_bvarPANELs_sample_nu`, aux_nu, aux_Sigma_c, aux_Sigma, prior) +} + sample_Sigma <- function(aux_Sigma_c_inv, aux_s, aux_nu, prior) { .Call(`_bvarPANELs_sample_Sigma`, aux_Sigma_c_inv, aux_s, aux_nu, prior) } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index c95a403..8fec022 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -84,6 +84,20 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// sample_nu +double sample_nu(const double& aux_nu, const arma::cube& aux_Sigma_c, const arma::mat& aux_Sigma, const Rcpp::List& prior); +RcppExport SEXP _bvarPANELs_sample_nu(SEXP aux_nuSEXP, SEXP aux_Sigma_cSEXP, SEXP aux_SigmaSEXP, SEXP priorSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const double& >::type aux_nu(aux_nuSEXP); + Rcpp::traits::input_parameter< const arma::cube& >::type aux_Sigma_c(aux_Sigma_cSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type aux_Sigma(aux_SigmaSEXP); + Rcpp::traits::input_parameter< const Rcpp::List& >::type prior(priorSEXP); + rcpp_result_gen = Rcpp::wrap(sample_nu(aux_nu, aux_Sigma_c, aux_Sigma, prior)); + return rcpp_result_gen; +END_RCPP +} // sample_Sigma arma::mat sample_Sigma(const arma::cube& aux_Sigma_c_inv, const double& aux_s, const double& aux_nu, const Rcpp::List& prior); RcppExport SEXP _bvarPANELs_sample_Sigma(SEXP aux_Sigma_c_invSEXP, SEXP aux_sSEXP, SEXP aux_nuSEXP, SEXP priorSEXP) { @@ -137,6 +151,7 @@ static const R_CallMethodDef CallEntries[] = { {"_bvarPANELs_sample_w", (DL_FUNC) &_bvarPANELs_sample_w, 2}, {"_bvarPANELs_sample_s", (DL_FUNC) &_bvarPANELs_sample_s, 5}, {"_bvarPANELs_log_kernel_nu", (DL_FUNC) &_bvarPANELs_log_kernel_nu, 7}, + {"_bvarPANELs_sample_nu", (DL_FUNC) &_bvarPANELs_sample_nu, 4}, {"_bvarPANELs_sample_Sigma", (DL_FUNC) &_bvarPANELs_sample_Sigma, 4}, {"_bvarPANELs_sample_AV", (DL_FUNC) &_bvarPANELs_sample_AV, 6}, {"_bvarPANELs_sample_A_c_Sigma_c", (DL_FUNC) &_bvarPANELs_sample_A_c_Sigma_c, 6}, diff --git a/src/sample_mniw.cpp b/src/sample_mniw.cpp index 18f81ac..4a2aee3 100644 --- a/src/sample_mniw.cpp +++ b/src/sample_mniw.cpp @@ -141,11 +141,45 @@ double log_kernel_nu ( - - - - - +// [[Rcpp:interface(cpp)]] +// [[Rcpp::export]] +double sample_nu ( + const double& aux_nu, // scalar + const arma::cube& aux_Sigma_c,// NxNxC + const arma::mat& aux_Sigma, // NxN + const Rcpp::List& prior +) { + + double prior_lambda = as(prior["lambda"]); + mat prior_M = as(prior["M"]); + int K = prior_M.n_rows; + int C = aux_Sigma_c.n_slices; + int N = aux_Sigma.n_rows; + + + double Cov_nu = 4 / C; + double Cov_nu_tmp = 0; + for (int n = 1; n < N + 1; n++) { + Cov_nu_tmp += R::psigamma( 0.5 * (aux_nu + 1 - n), 1); + } // END n loop + Cov_nu += 1 / Cov_nu_tmp; + + // Metropolis-Hastings + double aux_nu_star = RcppTN::rtn1( aux_nu, Cov_nu, N + 1, R_PosInf ); + double lk_nu_star = log_kernel_nu ( aux_nu_star, aux_Sigma_c, aux_Sigma, prior_lambda, C, N, K ); + double lk_nu_old = log_kernel_nu ( aux_nu, aux_Sigma_c, aux_Sigma, prior_lambda, C, N, K ); + double cgd_ratio = RcppTN::dtn1( aux_nu_star, aux_nu, Cov_nu, N + 1, R_PosInf ) / RcppTN::dtn1( aux_nu, aux_nu_star, Cov_nu, N + 1, R_PosInf ); + + double u = randu(); + double out = 0; + if (u < exp(lk_nu_star - lk_nu_old) * cgd_ratio) { + out = aux_nu_star; + } else { + out = aux_nu; + } // + + return out; +} // END sample_nu diff --git a/src/sample_mniw.h b/src/sample_mniw.h index bcb437c..b48f992 100644 --- a/src/sample_mniw.h +++ b/src/sample_mniw.h @@ -50,6 +50,14 @@ double log_kernel_nu ( ); +double sample_nu ( + const double& aux_nu, // scalar + const arma::cube& aux_Sigma_c,// NxNxC + const arma::mat& aux_Sigma, // NxN + const Rcpp::List& prior +); + + arma::mat sample_Sigma ( const arma::cube& aux_Sigma_c_inv, // NxNxC const double& aux_s, // scalar From 226776c540ba4e77cbe78f2375e88a5a54a03a28 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 5 May 2024 19:05:34 +1000 Subject: [PATCH 19/32] Include RcppProgress #2 --- DESCRIPTION | 7 ++++--- NAMESPACE | 1 + R/bvarPANELs-package.R | 1 + 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 9d4e43e..25da05f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -8,9 +8,10 @@ Author: c(person(given="Tomasz", family="Woźniak", email="wozniak.tom@pm.me", r Maintainer: Tomasz Woźniak License: GPL (>= 3) Imports: + R6, Rcpp (>= 1.0.12), - RcppTN, - R6 -LinkingTo: Rcpp, RcppArmadillo, RcppTN + RcppProgress, + RcppTN +LinkingTo: Rcpp, RcppArmadillo, RcppProgress, RcppTN Encoding: UTF-8 RoxygenNote: 7.3.1 diff --git a/NAMESPACE b/NAMESPACE index 0d3b7c0..2ad45d4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -5,6 +5,7 @@ export(specify_panel_data_matrices) export(specify_posterior_bvarPANEL) export(specify_prior_bvarPANEL) export(specify_starting_values_bvarPANEL) +import(RcppProgress) importFrom(R6,R6Class) importFrom(Rcpp,sourceCpp) importFrom(RcppTN,dtn) diff --git a/R/bvarPANELs-package.R b/R/bvarPANELs-package.R index 936837a..2b1335c 100644 --- a/R/bvarPANELs-package.R +++ b/R/bvarPANELs-package.R @@ -33,6 +33,7 @@ #' @importFrom Rcpp sourceCpp #' @importFrom R6 R6Class #' @importFrom RcppTN rtn dtn +#' @import RcppProgress #' @note This package is currently in active development. #' @author Tomasz Woźniak \email{wozniak.tom@pm.me} #' @keywords package models ts From f5e7723bb1a96d0186a980e75274ef999228ae7b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 5 May 2024 19:41:13 +1000 Subject: [PATCH 20/32] included ordinal in utils #2 --- src/utils.cpp | 24 ++++++++++++++++++++++++ src/utils.h | 15 +++++++++++++++ 2 files changed, 39 insertions(+) create mode 100644 src/utils.cpp create mode 100644 src/utils.h diff --git a/src/utils.cpp b/src/utils.cpp new file mode 100644 index 0000000..bf9d682 --- /dev/null +++ b/src/utils.cpp @@ -0,0 +1,24 @@ + +#include + +using namespace Rcpp; +using namespace arma; + + +// [[Rcpp::interfaces(cpp)]] +// [[Rcpp::export]] +std::string ordinal( + int n +) { + std::string suffix; + if (n % 10 == 1 && n % 100 != 11) { + suffix = "st"; + } else if (n % 10 == 2 && n % 100 != 12) { + suffix = "nd"; + } else if (n % 10 == 3 && n % 100 != 13) { + suffix = "rd"; + } else { + suffix = "th"; + } + return std::to_string(n) + suffix; +} // END ordinal diff --git a/src/utils.h b/src/utils.h new file mode 100644 index 0000000..da1b4ea --- /dev/null +++ b/src/utils.h @@ -0,0 +1,15 @@ +#ifndef _UTILS_H_ +#define _UTILS_H_ + +#include + +using namespace Rcpp; +using namespace arma; + + +std::string ordinal( + int n +); + + +#endif // _UTILS_H_ From 8e74575f47f9b63bd7937cbed31b10d35bb3c54c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 5 May 2024 19:42:11 +1000 Subject: [PATCH 21/32] a working bvarPANEL function structure #2 does no sampling ATM --- R/RcppExports.R | 8 ++ inst/include/bvarPANELs.h | 9 ++ inst/include/bvarPANELs_RcppExports.h | 51 +++++++++++ src/RcppExports.cpp | 71 +++++++++++++++ src/bvarPANEL.cpp | 120 ++++++++++++++++++++++++++ src/bvarPANEL.h | 19 ++++ 6 files changed, 278 insertions(+) create mode 100644 inst/include/bvarPANELs.h create mode 100644 inst/include/bvarPANELs_RcppExports.h create mode 100644 src/bvarPANEL.cpp create mode 100644 src/bvarPANEL.h diff --git a/R/RcppExports.R b/R/RcppExports.R index 8e6f27d..3b4e2fa 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,6 +1,10 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 +bvarPANEL <- function(S, prior, starting_values, thin = 100L, show_progress = TRUE) { + .Call(`_bvarPANELs_bvarPANEL`, S, prior, starting_values, thin, show_progress) +} + rmniw1 <- function(A, V, S, nu) { .Call(`_bvarPANELs_rmniw1`, A, V, S, nu) } @@ -37,3 +41,7 @@ sample_A_c_Sigma_c <- function(Y_c, X_c, aux_A, aux_V, aux_Sigma, aux_nu) { .Call(`_bvarPANELs_sample_A_c_Sigma_c`, Y_c, X_c, aux_A, aux_V, aux_Sigma, aux_nu) } +# Register entry points for exported C++ functions +methods::setLoadAction(function(ns) { + .Call(`_bvarPANELs_RcppExport_registerCCallable`) +}) diff --git a/inst/include/bvarPANELs.h b/inst/include/bvarPANELs.h new file mode 100644 index 0000000..e61a267 --- /dev/null +++ b/inst/include/bvarPANELs.h @@ -0,0 +1,9 @@ +// Generated by using Rcpp::compileAttributes() -> do not edit by hand +// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#ifndef RCPP_bvarPANELs_H_GEN_ +#define RCPP_bvarPANELs_H_GEN_ + +#include "bvarPANELs_RcppExports.h" + +#endif // RCPP_bvarPANELs_H_GEN_ diff --git a/inst/include/bvarPANELs_RcppExports.h b/inst/include/bvarPANELs_RcppExports.h new file mode 100644 index 0000000..7258d72 --- /dev/null +++ b/inst/include/bvarPANELs_RcppExports.h @@ -0,0 +1,51 @@ +// Generated by using Rcpp::compileAttributes() -> do not edit by hand +// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#ifndef RCPP_bvarPANELs_RCPPEXPORTS_H_GEN_ +#define RCPP_bvarPANELs_RCPPEXPORTS_H_GEN_ + +#include +#include + +namespace bvarPANELs { + + using namespace Rcpp; + + namespace { + void validateSignature(const char* sig) { + Rcpp::Function require = Rcpp::Environment::base_env()["require"]; + require("bvarPANELs", Rcpp::Named("quietly") = true); + typedef int(*Ptr_validate)(const char*); + static Ptr_validate p_validate = (Ptr_validate) + R_GetCCallable("bvarPANELs", "_bvarPANELs_RcppExport_validate"); + if (!p_validate(sig)) { + throw Rcpp::function_not_exported( + "C++ function with signature '" + std::string(sig) + "' not found in bvarPANELs"); + } + } + } + + inline std::string ordinal(int n) { + typedef SEXP(*Ptr_ordinal)(SEXP); + static Ptr_ordinal p_ordinal = NULL; + if (p_ordinal == NULL) { + validateSignature("std::string(*ordinal)(int)"); + p_ordinal = (Ptr_ordinal)R_GetCCallable("bvarPANELs", "_bvarPANELs_ordinal"); + } + RObject rcpp_result_gen; + { + RNGScope RCPP_rngScope_gen; + rcpp_result_gen = p_ordinal(Shield(Rcpp::wrap(n))); + } + if (rcpp_result_gen.inherits("interrupted-error")) + throw Rcpp::internal::InterruptedException(); + if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen)) + throw Rcpp::LongjumpException(rcpp_result_gen); + if (rcpp_result_gen.inherits("try-error")) + throw Rcpp::exception(Rcpp::as(rcpp_result_gen).c_str()); + return Rcpp::as(rcpp_result_gen); + } + +} + +#endif // RCPP_bvarPANELs_RCPPEXPORTS_H_GEN_ diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 8fec022..63ca6c8 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -1,8 +1,11 @@ // Generated by using Rcpp::compileAttributes() -> do not edit by hand // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 +#include "../inst/include/bvarPANELs.h" #include #include +#include +#include using namespace Rcpp; @@ -11,6 +14,21 @@ Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif +// bvarPANEL +Rcpp::List bvarPANEL(const int& S, const Rcpp::List& prior, const Rcpp::List& starting_values, const int thin, const bool show_progress); +RcppExport SEXP _bvarPANELs_bvarPANEL(SEXP SSEXP, SEXP priorSEXP, SEXP starting_valuesSEXP, SEXP thinSEXP, SEXP show_progressSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const int& >::type S(SSEXP); + Rcpp::traits::input_parameter< const Rcpp::List& >::type prior(priorSEXP); + Rcpp::traits::input_parameter< const Rcpp::List& >::type starting_values(starting_valuesSEXP); + Rcpp::traits::input_parameter< const int >::type thin(thinSEXP); + Rcpp::traits::input_parameter< const bool >::type show_progress(show_progressSEXP); + rcpp_result_gen = Rcpp::wrap(bvarPANEL(S, prior, starting_values, thin, show_progress)); + return rcpp_result_gen; +END_RCPP +} // rmniw1 arma::field rmniw1(const arma::mat& A, const arma::mat& V, const arma::mat& S, const double& nu); RcppExport SEXP _bvarPANELs_rmniw1(SEXP ASEXP, SEXP VSEXP, SEXP SSEXP, SEXP nuSEXP) { @@ -144,8 +162,59 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// ordinal +std::string ordinal(int n); +static SEXP _bvarPANELs_ordinal_try(SEXP nSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::traits::input_parameter< int >::type n(nSEXP); + rcpp_result_gen = Rcpp::wrap(ordinal(n)); + return rcpp_result_gen; +END_RCPP_RETURN_ERROR +} +RcppExport SEXP _bvarPANELs_ordinal(SEXP nSEXP) { + SEXP rcpp_result_gen; + { + Rcpp::RNGScope rcpp_rngScope_gen; + rcpp_result_gen = PROTECT(_bvarPANELs_ordinal_try(nSEXP)); + } + Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error"); + if (rcpp_isInterrupt_gen) { + UNPROTECT(1); + Rf_onintr(); + } + bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen); + if (rcpp_isLongjump_gen) { + Rcpp::internal::resumeJump(rcpp_result_gen); + } + Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error"); + if (rcpp_isError_gen) { + SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen); + UNPROTECT(1); + Rf_error("%s", CHAR(rcpp_msgSEXP_gen)); + } + UNPROTECT(1); + return rcpp_result_gen; +} + +// validate (ensure exported C++ functions exist before calling them) +static int _bvarPANELs_RcppExport_validate(const char* sig) { + static std::set signatures; + if (signatures.empty()) { + signatures.insert("std::string(*ordinal)(int)"); + } + return signatures.find(sig) != signatures.end(); +} + +// registerCCallable (register entry points for exported C++ functions) +RcppExport SEXP _bvarPANELs_RcppExport_registerCCallable() { + R_RegisterCCallable("bvarPANELs", "_bvarPANELs_ordinal", (DL_FUNC)_bvarPANELs_ordinal_try); + R_RegisterCCallable("bvarPANELs", "_bvarPANELs_RcppExport_validate", (DL_FUNC)_bvarPANELs_RcppExport_validate); + return R_NilValue; +} static const R_CallMethodDef CallEntries[] = { + {"_bvarPANELs_bvarPANEL", (DL_FUNC) &_bvarPANELs_bvarPANEL, 5}, {"_bvarPANELs_rmniw1", (DL_FUNC) &_bvarPANELs_rmniw1, 4}, {"_bvarPANELs_sample_m", (DL_FUNC) &_bvarPANELs_sample_m, 5}, {"_bvarPANELs_sample_w", (DL_FUNC) &_bvarPANELs_sample_w, 2}, @@ -155,6 +224,8 @@ static const R_CallMethodDef CallEntries[] = { {"_bvarPANELs_sample_Sigma", (DL_FUNC) &_bvarPANELs_sample_Sigma, 4}, {"_bvarPANELs_sample_AV", (DL_FUNC) &_bvarPANELs_sample_AV, 6}, {"_bvarPANELs_sample_A_c_Sigma_c", (DL_FUNC) &_bvarPANELs_sample_A_c_Sigma_c, 6}, + {"_bvarPANELs_ordinal", (DL_FUNC) &_bvarPANELs_ordinal, 1}, + {"_bvarPANELs_RcppExport_registerCCallable", (DL_FUNC) &_bvarPANELs_RcppExport_registerCCallable, 0}, {NULL, NULL, 0} }; diff --git a/src/bvarPANEL.cpp b/src/bvarPANEL.cpp new file mode 100644 index 0000000..28bf68d --- /dev/null +++ b/src/bvarPANEL.cpp @@ -0,0 +1,120 @@ +#include +#include "progress.hpp" + +#include "sample_mniw.h" +#include "utils.h" + +using namespace Rcpp; +using namespace arma; + + +// [[Rcpp:interface(cpp,r)]] +// [[Rcpp::export]] +Rcpp::List bvarPANEL( + const int& S, // No. of posterior draws + const Rcpp::List& prior, // a list of priors + const Rcpp::List& starting_values, + const int thin = 100, // introduce thinning + const bool show_progress = true +) { + // const Rcpp::List& Y, // a C-list of T_cxN elements + // const Rcpp::List& X, // a C-list of T_cxK elements + + // Progress bar setup + vec prog_rep_points = arma::round(arma::linspace(0, S, 50)); + + std::string oo = ""; + if ( thin != 1 ) { + oo = ordinal(thin) + " "; + } + + if (show_progress) { + Rcout << "**************************************************|" << endl; + Rcout << "bvarPANELs: Forecasting with Bayesian Hierarchical|" << endl; + Rcout << " Panel Vector Autoregressions |" << endl; + Rcout << "**************************************************|" << endl; + // Rcout << " Gibbs sampler for the SVAR-SV model |" << endl; + // Rcout << "**************************************************|" << endl; + Rcout << " Progress of the MCMC simulation for " << S << " draws" << endl; + Rcout << " Every " << oo << "draw is saved via MCMC thinning" << endl; + Rcout << " Press Esc to interrupt the computations" << endl; + Rcout << "**************************************************|" << endl; + } + Progress p(50, show_progress); + + cube aux_A_c = as(starting_values["A_c"]); + cube aux_Sigma_c = as(starting_values["Sigma_c"]); + mat aux_A = as(starting_values["A"]); + mat aux_V = as(starting_values["V"]); + mat aux_Sigma = as(starting_values["Sigma"]); + double aux_nu = as(starting_values["nu"]); + double aux_m = as(starting_values["m"]); + double aux_w = as(starting_values["w"]); + double aux_s = as(starting_values["s"]); + + const int C = aux_A_c.n_slices; + const int N = aux_A.n_cols; + const int K = aux_A.n_rows; + + field posterior_A_c(S); + field posterior_Sigma_c(S); + cube posterior_A(K, N, S); + cube posterior_V(K, K, S); + cube posterior_Sigma(N, N, S); + vec posterior_nu(S); + vec posterior_m(S); + vec posterior_w(S); + vec posterior_s(S); + + int ss = 0; + + for (int s=0; s + +using namespace Rcpp; +using namespace arma; + + +Rcpp::List bvarPANEL( + const int& S, // No. of posterior draws + const Rcpp::List& prior, // a list of priors + const Rcpp::List& starting_values, + const int thin = 100, // introduce thinning + const bool show_progress = true +); + + +#endif // _BVARPANEL_H_ From 16b5f4da274b7e552e88954296204799ab3d97fc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 5 May 2024 21:44:45 +1000 Subject: [PATCH 22/32] getting rid of utils #2 --- inst/include/bvarPANELs_RcppExports.h | 51 --------------------------- src/utils.cpp | 24 ------------- src/utils.h | 15 -------- 3 files changed, 90 deletions(-) delete mode 100644 inst/include/bvarPANELs_RcppExports.h delete mode 100644 src/utils.cpp delete mode 100644 src/utils.h diff --git a/inst/include/bvarPANELs_RcppExports.h b/inst/include/bvarPANELs_RcppExports.h deleted file mode 100644 index 7258d72..0000000 --- a/inst/include/bvarPANELs_RcppExports.h +++ /dev/null @@ -1,51 +0,0 @@ -// Generated by using Rcpp::compileAttributes() -> do not edit by hand -// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 - -#ifndef RCPP_bvarPANELs_RCPPEXPORTS_H_GEN_ -#define RCPP_bvarPANELs_RCPPEXPORTS_H_GEN_ - -#include -#include - -namespace bvarPANELs { - - using namespace Rcpp; - - namespace { - void validateSignature(const char* sig) { - Rcpp::Function require = Rcpp::Environment::base_env()["require"]; - require("bvarPANELs", Rcpp::Named("quietly") = true); - typedef int(*Ptr_validate)(const char*); - static Ptr_validate p_validate = (Ptr_validate) - R_GetCCallable("bvarPANELs", "_bvarPANELs_RcppExport_validate"); - if (!p_validate(sig)) { - throw Rcpp::function_not_exported( - "C++ function with signature '" + std::string(sig) + "' not found in bvarPANELs"); - } - } - } - - inline std::string ordinal(int n) { - typedef SEXP(*Ptr_ordinal)(SEXP); - static Ptr_ordinal p_ordinal = NULL; - if (p_ordinal == NULL) { - validateSignature("std::string(*ordinal)(int)"); - p_ordinal = (Ptr_ordinal)R_GetCCallable("bvarPANELs", "_bvarPANELs_ordinal"); - } - RObject rcpp_result_gen; - { - RNGScope RCPP_rngScope_gen; - rcpp_result_gen = p_ordinal(Shield(Rcpp::wrap(n))); - } - if (rcpp_result_gen.inherits("interrupted-error")) - throw Rcpp::internal::InterruptedException(); - if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen)) - throw Rcpp::LongjumpException(rcpp_result_gen); - if (rcpp_result_gen.inherits("try-error")) - throw Rcpp::exception(Rcpp::as(rcpp_result_gen).c_str()); - return Rcpp::as(rcpp_result_gen); - } - -} - -#endif // RCPP_bvarPANELs_RCPPEXPORTS_H_GEN_ diff --git a/src/utils.cpp b/src/utils.cpp deleted file mode 100644 index bf9d682..0000000 --- a/src/utils.cpp +++ /dev/null @@ -1,24 +0,0 @@ - -#include - -using namespace Rcpp; -using namespace arma; - - -// [[Rcpp::interfaces(cpp)]] -// [[Rcpp::export]] -std::string ordinal( - int n -) { - std::string suffix; - if (n % 10 == 1 && n % 100 != 11) { - suffix = "st"; - } else if (n % 10 == 2 && n % 100 != 12) { - suffix = "nd"; - } else if (n % 10 == 3 && n % 100 != 13) { - suffix = "rd"; - } else { - suffix = "th"; - } - return std::to_string(n) + suffix; -} // END ordinal diff --git a/src/utils.h b/src/utils.h deleted file mode 100644 index da1b4ea..0000000 --- a/src/utils.h +++ /dev/null @@ -1,15 +0,0 @@ -#ifndef _UTILS_H_ -#define _UTILS_H_ - -#include - -using namespace Rcpp; -using namespace arma; - - -std::string ordinal( - int n -); - - -#endif // _UTILS_H_ From ce04a7925b8002947a1b38113b0dc4f56fb2ed76 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 5 May 2024 21:45:42 +1000 Subject: [PATCH 23/32] bringing in bsvars #2 --- DESCRIPTION | 3 ++- R/bvarPANELs-package.R | 1 + src/bvarPANEL.cpp | 4 ++-- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 25da05f..666825d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -8,10 +8,11 @@ Author: c(person(given="Tomasz", family="Woźniak", email="wozniak.tom@pm.me", r Maintainer: Tomasz Woźniak License: GPL (>= 3) Imports: + bsvars, R6, Rcpp (>= 1.0.12), RcppProgress, RcppTN -LinkingTo: Rcpp, RcppArmadillo, RcppProgress, RcppTN +LinkingTo: bsvars, Rcpp, RcppArmadillo, RcppProgress, RcppTN Encoding: UTF-8 RoxygenNote: 7.3.1 diff --git a/R/bvarPANELs-package.R b/R/bvarPANELs-package.R index 2b1335c..5dd105a 100644 --- a/R/bvarPANELs-package.R +++ b/R/bvarPANELs-package.R @@ -30,6 +30,7 @@ #' @aliases bvarPANELs-package bvarPANELs #' @docType package #' @useDynLib bvarPANELs, .registration = TRUE +#' @importFrom bsvars estimate #' @importFrom Rcpp sourceCpp #' @importFrom R6 R6Class #' @importFrom RcppTN rtn dtn diff --git a/src/bvarPANEL.cpp b/src/bvarPANEL.cpp index 28bf68d..a186ffa 100644 --- a/src/bvarPANEL.cpp +++ b/src/bvarPANEL.cpp @@ -1,8 +1,8 @@ #include #include "progress.hpp" +#include "bsvars.h" #include "sample_mniw.h" -#include "utils.h" using namespace Rcpp; using namespace arma; @@ -25,7 +25,7 @@ Rcpp::List bvarPANEL( std::string oo = ""; if ( thin != 1 ) { - oo = ordinal(thin) + " "; + oo = bsvars::ordinal(thin) + " "; } if (show_progress) { From e01e78857f3a2ce14f3d59b539877ae4809ee139 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 5 May 2024 21:46:27 +1000 Subject: [PATCH 24/32] necessary corrections to specify #2 to make a basic example with data work --- R/specify_bvarpanel.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/specify_bvarpanel.R b/R/specify_bvarpanel.R index 5c5e36d..389ce6b 100644 --- a/R/specify_bvarpanel.R +++ b/R/specify_bvarpanel.R @@ -205,7 +205,7 @@ specify_starting_values_bvarPANEL = R6::R6Class( stopifnot("Argument p must be a positive integer number." = p > 0 & p %% 1 == 0) K = N * p + 1 - self$A_c = matrix(stats::rnorm(C * K * N, sd = 0.001), c(K, N, C)) + self$A_c = array(stats::rnorm(C * K * N, sd = 0.001), c(K, N, C)) self$Sigma_c = stats::rWishart(C, N + 1, diag(N)) self$A = matrix(stats::rnorm(K * N, sd = 0.001), K, N) + diag(K)[,1:N] self$V = stats::rWishart(1, K + 1, diag(K))[,,1] @@ -306,7 +306,7 @@ specify_panel_data_matrices = R6::R6Class( } else { stopifnot("Argument data has to be a list of matrices." = is.list(data) & all(simplify2array(lapply(data, function(x){is.matrix(x) & is.numeric(x)})))) stopifnot("Argument data has to contain matrices with the same number of columns." = length(unique(simplify2array(lapply(data, ncol)))) == 1) - stopifnot("Argument data cannot include missing values." = all(simplify2array(lapply(dd, function(x){!any(is.na(x))})))) + stopifnot("Argument data cannot include missing values." = all(simplify2array(lapply(data, function(x){!any(is.na(x))})))) } stopifnot("Argument p must be a positive integer number." = p > 0 & p %% 1 == 0) @@ -321,7 +321,7 @@ specify_panel_data_matrices = R6::R6Class( X = cbind(X, data[[c]][(p + 1):TT - i,]) } X = cbind(X, rep(1, T_c)) - self$X = X + self$X[[c]] = X } # END c loop names(self$Y) = names(self$X) = names(data) }, # END initialize From 7672541a92d2c05780095781c4ca030e337d7eb4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 5 May 2024 21:47:06 +1000 Subject: [PATCH 25/32] creating estimate methods #2 --- NAMESPACE | 3 ++ R/RcppExports.R | 4 -- R/estimate.bvarPANEL.R | 87 ++++++++++++++++++++++++++++++ inst/include/bvarPANELs.h | 9 ---- man/estimate.BVARPANEL.Rd | 49 +++++++++++++++++ man/estimate.PosteriorBVARPANEL.Rd | 51 ++++++++++++++++++ src/RcppExports.cpp | 55 ------------------- 7 files changed, 190 insertions(+), 68 deletions(-) create mode 100644 R/estimate.bvarPANEL.R delete mode 100644 inst/include/bvarPANELs.h create mode 100644 man/estimate.BVARPANEL.Rd create mode 100644 man/estimate.PosteriorBVARPANEL.Rd diff --git a/NAMESPACE b/NAMESPACE index 2ad45d4..7902c2d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,7 @@ # Generated by roxygen2: do not edit by hand +S3method(estimate,BVARPANEL) +S3method(estimate,PosteriorBVARPANEL) export(specify_bvarPANEL) export(specify_panel_data_matrices) export(specify_posterior_bvarPANEL) @@ -10,4 +12,5 @@ importFrom(R6,R6Class) importFrom(Rcpp,sourceCpp) importFrom(RcppTN,dtn) importFrom(RcppTN,rtn) +importFrom(bsvars,estimate) useDynLib(bvarPANELs, .registration = TRUE) diff --git a/R/RcppExports.R b/R/RcppExports.R index 3b4e2fa..ae1ee32 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -41,7 +41,3 @@ sample_A_c_Sigma_c <- function(Y_c, X_c, aux_A, aux_V, aux_Sigma, aux_nu) { .Call(`_bvarPANELs_sample_A_c_Sigma_c`, Y_c, X_c, aux_A, aux_V, aux_Sigma, aux_nu) } -# Register entry points for exported C++ functions -methods::setLoadAction(function(ns) { - .Call(`_bvarPANELs_RcppExport_registerCCallable`) -}) diff --git a/R/estimate.bvarPANEL.R b/R/estimate.bvarPANEL.R new file mode 100644 index 0000000..f53a27b --- /dev/null +++ b/R/estimate.bvarPANEL.R @@ -0,0 +1,87 @@ + +#' @title Bayesian estimation of a Bayesian Hierarchical Panel Vector +#' Autoregression using Gibbs sampler +#' +#' @description Estimates the Bayesian Hierarchical Panel VAR using the Gibbs +#' sampler proposed by Woźniak (2024). +#' +#' @details +#' The homoskedastic SVAR model is given by the reduced form equation: +#' \deqn{Y_c = A_cX_c + E_c} +#' where \eqn{Y_c} is an \code{T_c x N} matrix of dependent variables for +#' country \code{c}, \eqn{X_c} is a \code{T_c x K} matrix of explanatory +#' variables, \eqn{E_c} is an \code{T_c x N} matrix of error terms, and +#' \eqn{A_c} is an \code{NxK} matrix of country-specific autoregressive slope +#' coefficients and parameters on deterministic terms in \eqn{X_c}. +#' +#' @param specification an object of class BVARPANEL generated using the +#' \code{specify_bvarPANEL$new()} function. +#' @param S a positive integer, the number of posterior draws to be generated +#' @param thin a positive integer, specifying the frequency of MCMC output thinning +#' @param show_progress a logical value, if \code{TRUE} the estimation progress +#' bar is visible +#' +#' @return An object of class PosteriorBVARPANEL containing the Bayesian +#' estimation output and containing two elements: +#' +#' \code{posterior} a list with a collection of \code{S} draws from the +#' posterior distribution generated via Gibbs sampler. +#' \code{last_draw} an object of class BVARPANEL with the last draw of the +#' current MCMC run as the starting value to be passed to the continuation of +#' the MCMC estimation using the \code{estimate()} method. +#' +#' @seealso \code{\link{specify_bvarPANEL}}, \code{\link{specify_posterior_bvarPANEL}} +#' +#' @author Tomasz Woźniak \email{wozniak.tom@pm.me} +#' +#' @method estimate BVARPANEL +#' +#' @export +estimate.BVARPANEL <- function( + specification, + S, + thin = 10L, + show_progress = TRUE +) { + + # get the inputs to estimation + prior = specification$prior$get_prior() + starting_values = specification$starting_values$get_starting_values() + # data_matrices = specification$data_matrices$get_data_matrices() + + # estimation + qqq = .Call(`_bvarPANELs_bvarPANEL`, S, prior, starting_values, thin, show_progress) + + specification$starting_values$set_starting_values(qqq$last_draw) + output = specify_posterior_bvarPANEL$new(specification, qqq$posterior) + + return(output) +} # END estimate.BVARPANEL + + + +#' @inherit estimate.BVARPANEL +#' +#' @method estimate PosteriorBVARPANEL +#' +#' @param specification an object of class PosteriorBVARPANEL generated using +#' the \code{estimate.BVARPANEL()} function. This setup facilitates the +#' continuation of the MCMC sampling starting from the last draw of the previous +#' run. +#' +#' @export +estimate.PosteriorBVARPANEL <- function(specification, S, thin = 10, show_progress = TRUE) { + + # get the inputs to estimation + prior = specification$last_draw$prior$get_prior() + starting_values = specification$last_draw$starting_values$get_starting_values() + # data_matrices = specification$last_draw$data_matrices$get_data_matrices() + + # estimation + qqq = .Call(`_bvarPANELs_bvarPANEL`, S, prior, starting_values, thin, show_progress) + + specification$starting_values$set_starting_values(qqq$last_draw) + output = specify_posterior_bvarPANEL$new(specification, qqq$posterior) + + return(output) +} # END estimate.PosteriorBSVAR diff --git a/inst/include/bvarPANELs.h b/inst/include/bvarPANELs.h deleted file mode 100644 index e61a267..0000000 --- a/inst/include/bvarPANELs.h +++ /dev/null @@ -1,9 +0,0 @@ -// Generated by using Rcpp::compileAttributes() -> do not edit by hand -// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 - -#ifndef RCPP_bvarPANELs_H_GEN_ -#define RCPP_bvarPANELs_H_GEN_ - -#include "bvarPANELs_RcppExports.h" - -#endif // RCPP_bvarPANELs_H_GEN_ diff --git a/man/estimate.BVARPANEL.Rd b/man/estimate.BVARPANEL.Rd new file mode 100644 index 0000000..111824f --- /dev/null +++ b/man/estimate.BVARPANEL.Rd @@ -0,0 +1,49 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/estimate.bvarPANEL.R +\name{estimate.BVARPANEL} +\alias{estimate.BVARPANEL} +\title{Bayesian estimation of a Bayesian Hierarchical Panel Vector +Autoregression using Gibbs sampler} +\usage{ +\method{estimate}{BVARPANEL}(specification, S, thin = 10L, show_progress = TRUE) +} +\arguments{ +\item{specification}{an object of class BVARPANEL generated using the +\code{specify_bvarPANEL$new()} function.} + +\item{S}{a positive integer, the number of posterior draws to be generated} + +\item{thin}{a positive integer, specifying the frequency of MCMC output thinning} + +\item{show_progress}{a logical value, if \code{TRUE} the estimation progress +bar is visible} +} +\value{ +An object of class PosteriorBVARPANEL containing the Bayesian +estimation output and containing two elements: + + \code{posterior} a list with a collection of \code{S} draws from the + posterior distribution generated via Gibbs sampler. +\code{last_draw} an object of class BVARPANEL with the last draw of the +current MCMC run as the starting value to be passed to the continuation of +the MCMC estimation using the \code{estimate()} method. +} +\description{ +Estimates the Bayesian Hierarchical Panel VAR using the Gibbs +sampler proposed by Woźniak (2024). +} +\details{ +The homoskedastic SVAR model is given by the reduced form equation: +\deqn{Y_c = A_cX_c + E_c} +where \eqn{Y_c} is an \code{T_c x N} matrix of dependent variables for +country \code{c}, \eqn{X_c} is a \code{T_c x K} matrix of explanatory +variables, \eqn{E_c} is an \code{T_c x N} matrix of error terms, and +\eqn{A_c} is an \code{NxK} matrix of country-specific autoregressive slope +coefficients and parameters on deterministic terms in \eqn{X_c}. +} +\seealso{ +\code{\link{specify_bvarPANEL}}, \code{\link{specify_posterior_bvarPANEL}} +} +\author{ +Tomasz Woźniak \email{wozniak.tom@pm.me} +} diff --git a/man/estimate.PosteriorBVARPANEL.Rd b/man/estimate.PosteriorBVARPANEL.Rd new file mode 100644 index 0000000..7bad6df --- /dev/null +++ b/man/estimate.PosteriorBVARPANEL.Rd @@ -0,0 +1,51 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/estimate.bvarPANEL.R +\name{estimate.PosteriorBVARPANEL} +\alias{estimate.PosteriorBVARPANEL} +\title{Bayesian estimation of a Bayesian Hierarchical Panel Vector +Autoregression using Gibbs sampler} +\usage{ +\method{estimate}{PosteriorBVARPANEL}(specification, S, thin = 10, show_progress = TRUE) +} +\arguments{ +\item{specification}{an object of class PosteriorBVARPANEL generated using +the \code{estimate.BVARPANEL()} function. This setup facilitates the +continuation of the MCMC sampling starting from the last draw of the previous +run.} + +\item{S}{a positive integer, the number of posterior draws to be generated} + +\item{thin}{a positive integer, specifying the frequency of MCMC output thinning} + +\item{show_progress}{a logical value, if \code{TRUE} the estimation progress +bar is visible} +} +\value{ +An object of class PosteriorBVARPANEL containing the Bayesian +estimation output and containing two elements: + + \code{posterior} a list with a collection of \code{S} draws from the + posterior distribution generated via Gibbs sampler. +\code{last_draw} an object of class BVARPANEL with the last draw of the +current MCMC run as the starting value to be passed to the continuation of +the MCMC estimation using the \code{estimate()} method. +} +\description{ +Estimates the Bayesian Hierarchical Panel VAR using the Gibbs +sampler proposed by Woźniak (2024). +} +\details{ +The homoskedastic SVAR model is given by the reduced form equation: +\deqn{Y_c = A_cX_c + E_c} +where \eqn{Y_c} is an \code{T_c x N} matrix of dependent variables for +country \code{c}, \eqn{X_c} is a \code{T_c x K} matrix of explanatory +variables, \eqn{E_c} is an \code{T_c x N} matrix of error terms, and +\eqn{A_c} is an \code{NxK} matrix of country-specific autoregressive slope +coefficients and parameters on deterministic terms in \eqn{X_c}. +} +\seealso{ +\code{\link{specify_bvarPANEL}}, \code{\link{specify_posterior_bvarPANEL}} +} +\author{ +Tomasz Woźniak \email{wozniak.tom@pm.me} +} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 63ca6c8..f11ecb8 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -1,11 +1,8 @@ // Generated by using Rcpp::compileAttributes() -> do not edit by hand // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 -#include "../inst/include/bvarPANELs.h" #include #include -#include -#include using namespace Rcpp; @@ -162,56 +159,6 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// ordinal -std::string ordinal(int n); -static SEXP _bvarPANELs_ordinal_try(SEXP nSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::traits::input_parameter< int >::type n(nSEXP); - rcpp_result_gen = Rcpp::wrap(ordinal(n)); - return rcpp_result_gen; -END_RCPP_RETURN_ERROR -} -RcppExport SEXP _bvarPANELs_ordinal(SEXP nSEXP) { - SEXP rcpp_result_gen; - { - Rcpp::RNGScope rcpp_rngScope_gen; - rcpp_result_gen = PROTECT(_bvarPANELs_ordinal_try(nSEXP)); - } - Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error"); - if (rcpp_isInterrupt_gen) { - UNPROTECT(1); - Rf_onintr(); - } - bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen); - if (rcpp_isLongjump_gen) { - Rcpp::internal::resumeJump(rcpp_result_gen); - } - Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error"); - if (rcpp_isError_gen) { - SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen); - UNPROTECT(1); - Rf_error("%s", CHAR(rcpp_msgSEXP_gen)); - } - UNPROTECT(1); - return rcpp_result_gen; -} - -// validate (ensure exported C++ functions exist before calling them) -static int _bvarPANELs_RcppExport_validate(const char* sig) { - static std::set signatures; - if (signatures.empty()) { - signatures.insert("std::string(*ordinal)(int)"); - } - return signatures.find(sig) != signatures.end(); -} - -// registerCCallable (register entry points for exported C++ functions) -RcppExport SEXP _bvarPANELs_RcppExport_registerCCallable() { - R_RegisterCCallable("bvarPANELs", "_bvarPANELs_ordinal", (DL_FUNC)_bvarPANELs_ordinal_try); - R_RegisterCCallable("bvarPANELs", "_bvarPANELs_RcppExport_validate", (DL_FUNC)_bvarPANELs_RcppExport_validate); - return R_NilValue; -} static const R_CallMethodDef CallEntries[] = { {"_bvarPANELs_bvarPANEL", (DL_FUNC) &_bvarPANELs_bvarPANEL, 5}, @@ -224,8 +171,6 @@ static const R_CallMethodDef CallEntries[] = { {"_bvarPANELs_sample_Sigma", (DL_FUNC) &_bvarPANELs_sample_Sigma, 4}, {"_bvarPANELs_sample_AV", (DL_FUNC) &_bvarPANELs_sample_AV, 6}, {"_bvarPANELs_sample_A_c_Sigma_c", (DL_FUNC) &_bvarPANELs_sample_A_c_Sigma_c, 6}, - {"_bvarPANELs_ordinal", (DL_FUNC) &_bvarPANELs_ordinal, 1}, - {"_bvarPANELs_RcppExport_registerCCallable", (DL_FUNC) &_bvarPANELs_RcppExport_registerCCallable, 0}, {NULL, NULL, 0} }; From ebdf618aa84b716d5e989a48cf86e0493f8a0d10 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sun, 5 May 2024 23:26:47 +1000 Subject: [PATCH 26/32] Gibbs sampler is working nicely! #2 --- R/RcppExports.R | 4 ++-- R/estimate.bvarPANEL.R | 12 +++++------ R/specify_bvarpanel.R | 2 +- src/RcppExports.cpp | 10 +++++---- src/bvarPANEL.cpp | 47 +++++++++++++++++++++++++++++++++++++++--- src/bvarPANEL.h | 2 ++ src/sample_mniw.cpp | 10 ++++++++- 7 files changed, 70 insertions(+), 17 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index ae1ee32..9649d98 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,8 +1,8 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 -bvarPANEL <- function(S, prior, starting_values, thin = 100L, show_progress = TRUE) { - .Call(`_bvarPANELs_bvarPANEL`, S, prior, starting_values, thin, show_progress) +bvarPANEL <- function(S, Y, X, prior, starting_values, thin = 100L, show_progress = TRUE) { + .Call(`_bvarPANELs_bvarPANEL`, S, Y, X, prior, starting_values, thin, show_progress) } rmniw1 <- function(A, V, S, nu) { diff --git a/R/estimate.bvarPANEL.R b/R/estimate.bvarPANEL.R index f53a27b..6f58991 100644 --- a/R/estimate.bvarPANEL.R +++ b/R/estimate.bvarPANEL.R @@ -47,10 +47,10 @@ estimate.BVARPANEL <- function( # get the inputs to estimation prior = specification$prior$get_prior() starting_values = specification$starting_values$get_starting_values() - # data_matrices = specification$data_matrices$get_data_matrices() + data_matrices = specification$data_matrices$get_data_matrices() # estimation - qqq = .Call(`_bvarPANELs_bvarPANEL`, S, prior, starting_values, thin, show_progress) + qqq = .Call(`_bvarPANELs_bvarPANEL`, S, data_matrices$Y, data_matrices$X, prior, starting_values, thin, show_progress) specification$starting_values$set_starting_values(qqq$last_draw) output = specify_posterior_bvarPANEL$new(specification, qqq$posterior) @@ -75,13 +75,13 @@ estimate.PosteriorBVARPANEL <- function(specification, S, thin = 10, show_progre # get the inputs to estimation prior = specification$last_draw$prior$get_prior() starting_values = specification$last_draw$starting_values$get_starting_values() - # data_matrices = specification$last_draw$data_matrices$get_data_matrices() + data_matrices = specification$last_draw$data_matrices$get_data_matrices() # estimation - qqq = .Call(`_bvarPANELs_bvarPANEL`, S, prior, starting_values, thin, show_progress) + qqq = .Call(`_bvarPANELs_bvarPANEL`, S, data_matrices$Y, data_matrices$X, prior, starting_values, thin, show_progress) - specification$starting_values$set_starting_values(qqq$last_draw) - output = specify_posterior_bvarPANEL$new(specification, qqq$posterior) + specification$last_draw$starting_values$set_starting_values(qqq$last_draw) + output = specify_posterior_bvarPANEL$new(specification$last_draw, qqq$posterior) return(output) } # END estimate.PosteriorBSVAR diff --git a/R/specify_bvarpanel.R b/R/specify_bvarpanel.R index 389ce6b..deb6dd7 100644 --- a/R/specify_bvarpanel.R +++ b/R/specify_bvarpanel.R @@ -498,7 +498,7 @@ specify_posterior_bvarPANEL = R6::R6Class( initialize = function(specification_bvarPANEL, posterior_bvarPANEL) { stopifnot("Argument specification_bvarPANEL must be of class BVARPANEL." = any(class(specification_bvarPANEL) == "BVARPANEL")) - stopifnot("Argument posterior_bvarPANEL must must contain MCMC output." = is.list(posterior_bvarPANEL) & is.array(posterior_bvarPANEL$A) & is.array(posterior_bvarPANEL$Sigma) & is.vector(posterior_bvarPANEL$w)) + stopifnot("Argument posterior_bvarPANEL must must contain MCMC output." = is.list(posterior_bvarPANEL) & is.array(posterior_bvarPANEL$A) & is.array(posterior_bvarPANEL$Sigma) & is.array(posterior_bvarPANEL$V)) self$last_draw = specification_bvarPANEL self$posterior = posterior_bvarPANEL diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index f11ecb8..d09c9cc 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -12,17 +12,19 @@ Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif // bvarPANEL -Rcpp::List bvarPANEL(const int& S, const Rcpp::List& prior, const Rcpp::List& starting_values, const int thin, const bool show_progress); -RcppExport SEXP _bvarPANELs_bvarPANEL(SEXP SSEXP, SEXP priorSEXP, SEXP starting_valuesSEXP, SEXP thinSEXP, SEXP show_progressSEXP) { +Rcpp::List bvarPANEL(const int& S, const Rcpp::List& Y, const Rcpp::List& X, const Rcpp::List& prior, const Rcpp::List& starting_values, const int thin, const bool show_progress); +RcppExport SEXP _bvarPANELs_bvarPANEL(SEXP SSEXP, SEXP YSEXP, SEXP XSEXP, SEXP priorSEXP, SEXP starting_valuesSEXP, SEXP thinSEXP, SEXP show_progressSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const int& >::type S(SSEXP); + Rcpp::traits::input_parameter< const Rcpp::List& >::type Y(YSEXP); + Rcpp::traits::input_parameter< const Rcpp::List& >::type X(XSEXP); Rcpp::traits::input_parameter< const Rcpp::List& >::type prior(priorSEXP); Rcpp::traits::input_parameter< const Rcpp::List& >::type starting_values(starting_valuesSEXP); Rcpp::traits::input_parameter< const int >::type thin(thinSEXP); Rcpp::traits::input_parameter< const bool >::type show_progress(show_progressSEXP); - rcpp_result_gen = Rcpp::wrap(bvarPANEL(S, prior, starting_values, thin, show_progress)); + rcpp_result_gen = Rcpp::wrap(bvarPANEL(S, Y, X, prior, starting_values, thin, show_progress)); return rcpp_result_gen; END_RCPP } @@ -161,7 +163,7 @@ END_RCPP } static const R_CallMethodDef CallEntries[] = { - {"_bvarPANELs_bvarPANEL", (DL_FUNC) &_bvarPANELs_bvarPANEL, 5}, + {"_bvarPANELs_bvarPANEL", (DL_FUNC) &_bvarPANELs_bvarPANEL, 7}, {"_bvarPANELs_rmniw1", (DL_FUNC) &_bvarPANELs_rmniw1, 4}, {"_bvarPANELs_sample_m", (DL_FUNC) &_bvarPANELs_sample_m, 5}, {"_bvarPANELs_sample_w", (DL_FUNC) &_bvarPANELs_sample_w, 2}, diff --git a/src/bvarPANEL.cpp b/src/bvarPANEL.cpp index a186ffa..52de5e3 100644 --- a/src/bvarPANEL.cpp +++ b/src/bvarPANEL.cpp @@ -12,13 +12,13 @@ using namespace arma; // [[Rcpp::export]] Rcpp::List bvarPANEL( const int& S, // No. of posterior draws + const Rcpp::List& Y, // a C-list of T_cxN elements + const Rcpp::List& X, // a C-list of T_cxK elements const Rcpp::List& prior, // a list of priors const Rcpp::List& starting_values, const int thin = 100, // introduce thinning const bool show_progress = true ) { - // const Rcpp::List& Y, // a C-list of T_cxN elements - // const Rcpp::List& X, // a C-list of T_cxK elements // Progress bar setup vec prog_rep_points = arma::round(arma::linspace(0, S, 50)); @@ -66,17 +66,58 @@ Rcpp::List bvarPANEL( vec posterior_w(S); vec posterior_s(S); + field y(C); + field x(C); + cube aux_Sigma_c_inv(N, N, C); + for (int c=0; c(Y[c]); + x(c) = as(X[c]); + aux_Sigma_c_inv.slice(c) = inv_sympd( aux_Sigma_c.slice(c) ); + } // END c loop + int ss = 0; for (int s=0; s tmp_AV = sample_AV( aux_A_c, aux_Sigma_c_inv, aux_s, aux_m, aux_w, prior ); + aux_A = tmp_AV(0); + aux_V = tmp_AV(1); + // sample aux_A_c, aux_Sigma_c + // Rcout << " sample A_c Sigma_c" << endl; + // Rcout << " aux_nu: " << aux_nu << endl; + for (int c=0; c tmp_A_c_Sigma_c = sample_A_c_Sigma_c( y(c), x(c), aux_A, aux_V, aux_Sigma, aux_nu ); + aux_A_c.slice(c) = tmp_A_c_Sigma_c(0); + aux_Sigma_c.slice(c) = tmp_A_c_Sigma_c(1); + aux_Sigma_c_inv.slice(c) = inv_sympd( aux_Sigma_c.slice(c) ); + } // END c loop if (s % thin == 0) { posterior_A_c(ss) = aux_A_c; diff --git a/src/bvarPANEL.h b/src/bvarPANEL.h index 79b690a..297e1c7 100644 --- a/src/bvarPANEL.h +++ b/src/bvarPANEL.h @@ -9,6 +9,8 @@ using namespace arma; Rcpp::List bvarPANEL( const int& S, // No. of posterior draws + const Rcpp::List& Y, // a C-list of T_cxN elements + const Rcpp::List& X, // a C-list of T_cxK elements const Rcpp::List& prior, // a list of priors const Rcpp::List& starting_values, const int thin = 100, // introduce thinning diff --git a/src/sample_mniw.cpp b/src/sample_mniw.cpp index 4a2aee3..4eee3c4 100644 --- a/src/sample_mniw.cpp +++ b/src/sample_mniw.cpp @@ -13,7 +13,13 @@ arma::field rmniw1( const arma::mat& S, // NxN const double& nu // scalar ) { - mat Sigma = iwishrnd(S, nu); + // Rcout << " A: " << A << std::endl; + // Rcout << " V: " << V.is_sympd() << std::endl; + // Rcout << " S: " << S.is_sympd() << std::endl; + // Rcout << " nu: " << nu << std::endl; + + mat SS = 0.5 * (S + S.t()); + mat Sigma = iwishrnd(SS, nu); mat X_tmp = mat(size(A), fill::randn); mat X = A + chol(V).t() * X_tmp * chol(Sigma); @@ -272,6 +278,8 @@ arma::field sample_A_c_Sigma_c ( mat Sigma_bar = aux_Sigma + Y_c.t() * Y_c + aux_A.t() * aux_V_inv * aux_A - A_bar.t() * V_bar_inv * A_bar; double nu_bar = T_c + aux_nu; + // Rcout << " nu_bar: " << nu_bar << std::endl; + // Rcout << " Sigma_bar: " << Sigma_bar.is_sympd() << std::endl; arma::field aux_A_c_Sigma_c = rmniw1( A_bar, V_bar, Sigma_bar, nu_bar ); return aux_A_c_Sigma_c; } // END sample_A_c_Sigma_c From 2d514cdbd4a419096fcf40cc45637ad2bef3dc8e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sat, 11 May 2024 22:48:42 +1000 Subject: [PATCH 27/32] a new specify help doc --- man/specify_bvarPANEL.Rd | 2 +- man/specify_panel_data_matrices.Rd | 2 +- man/specify_posterior_bvarPANEL.Rd | 2 +- man/specify_prior_bvarPANEL.Rd | 2 +- man/specify_starting_values_bvarPANEL.Rd | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/man/specify_bvarPANEL.Rd b/man/specify_bvarPANEL.Rd index 4d77395..8d7e888 100644 --- a/man/specify_bvarPANEL.Rd +++ b/man/specify_bvarPANEL.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/specify_bvarPANEL.R +% Please edit documentation in R/specify_bvarpanel.R \name{specify_bvarPANEL} \alias{specify_bvarPANEL} \title{R6 Class representing the specification of the BVARPANEL model} diff --git a/man/specify_panel_data_matrices.Rd b/man/specify_panel_data_matrices.Rd index c3242ae..284a246 100644 --- a/man/specify_panel_data_matrices.Rd +++ b/man/specify_panel_data_matrices.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/specify_bvarPANEL.R +% Please edit documentation in R/specify_bvarpanel.R \name{specify_panel_data_matrices} \alias{specify_panel_data_matrices} \title{R6 Class Representing DataMatricesBVARPANEL} diff --git a/man/specify_posterior_bvarPANEL.Rd b/man/specify_posterior_bvarPANEL.Rd index f361cd2..4a2ebe2 100644 --- a/man/specify_posterior_bvarPANEL.Rd +++ b/man/specify_posterior_bvarPANEL.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/specify_bvarPANEL.R +% Please edit documentation in R/specify_bvarpanel.R \name{specify_posterior_bvarPANEL} \alias{specify_posterior_bvarPANEL} \title{R6 Class Representing PosteriorBVARPANEL} diff --git a/man/specify_prior_bvarPANEL.Rd b/man/specify_prior_bvarPANEL.Rd index ac0dc8c..cb03743 100644 --- a/man/specify_prior_bvarPANEL.Rd +++ b/man/specify_prior_bvarPANEL.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/specify_bvarPANEL.R +% Please edit documentation in R/specify_bvarpanel.R \name{specify_prior_bvarPANEL} \alias{specify_prior_bvarPANEL} \title{R6 Class Representing PriorBVARPANEL} diff --git a/man/specify_starting_values_bvarPANEL.Rd b/man/specify_starting_values_bvarPANEL.Rd index 57a868d..8c821d9 100644 --- a/man/specify_starting_values_bvarPANEL.Rd +++ b/man/specify_starting_values_bvarPANEL.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/specify_bvarPANEL.R +% Please edit documentation in R/specify_bvarpanel.R \name{specify_starting_values_bvarPANEL} \alias{specify_starting_values_bvarPANEL} \title{R6 Class Representing StartingValuesBVARPANEL} From 1cd346c0f9195dc4847e55a4d10121f09adf71f3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sat, 11 May 2024 23:30:03 +1000 Subject: [PATCH 28/32] bsvars in Depends to use its generics #2 #3 --- DESCRIPTION | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index ea76ef7..171bf5e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -7,8 +7,10 @@ Date: 2024-04-27 Author: c(person(given="Tomasz", family="Woźniak", email="wozniak.tom@pm.me", role = c("aut", "cre"), comment = c(ORCID = "0000-0003-2212-2378")), person(given="Miguel", family="Sanchez-Martinez", role = "ctb"), person(given="International Labour Organization", role = "cph")) Maintainer: Tomasz Woźniak License: GPL (>= 3) +Depends: + R (>= 2.10), + bsvars Imports: - bsvars, R6, Rcpp (>= 1.0.12), RcppProgress, @@ -17,5 +19,3 @@ LinkingTo: bsvars, Rcpp, RcppArmadillo, RcppProgress, RcppTN LazyData: true Encoding: UTF-8 RoxygenNote: 7.3.1 -Depends: - R (>= 2.10) From 93bae8b66f8cc515199b56f5cbc594df2bf54b13 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sat, 11 May 2024 23:30:31 +1000 Subject: [PATCH 29/32] data allows full examples #3 --- R/specify_bvarpanel.R | 88 +++++++++++++++++------------- man/specify_bvarPANEL.Rd | 60 +++++++++----------- man/specify_panel_data_matrices.Rd | 26 +++++++++ man/specify_posterior_bvarPANEL.Rd | 69 +++++++++++++---------- 4 files changed, 140 insertions(+), 103 deletions(-) diff --git a/R/specify_bvarpanel.R b/R/specify_bvarpanel.R index deb6dd7..c3ae482 100644 --- a/R/specify_bvarpanel.R +++ b/R/specify_bvarpanel.R @@ -271,8 +271,6 @@ specify_starting_values_bvarPANEL = R6::R6Class( - - #' R6 Class Representing DataMatricesBVARPANEL #' #' @description @@ -280,6 +278,11 @@ specify_starting_values_bvarPANEL = R6::R6Class( #' variables, \eqn{\mathbf{Y}_c}, and regressors, \eqn{\mathbf{X}_c}, for the #' Bayesian Panel VAR model for all countries \eqn{c = 1, ..., C}. #' +#' @examples +#' data(ilo_cubic_panel) +#' YX = specify_panel_data_matrices$new(data = ilo_cubic_panel, p = 4) +#' length(YX$Y); names(YX$Y) +#' #' @export specify_panel_data_matrices = R6::R6Class( "DataMatricesBVARPANEL", @@ -328,6 +331,12 @@ specify_panel_data_matrices = R6::R6Class( #' @description #' Returns the data matrices DataMatricesBVARPANEL as a \code{list}. + #' + #' @examples + #' data(ilo_cubic_panel) + #' YX = specify_panel_data_matrices$new(ilo_cubic_panel) + #' YX$get_data_matrices() + #' get_data_matrices = function() { list( Y = self$Y, @@ -350,13 +359,11 @@ specify_panel_data_matrices = R6::R6Class( #' Vector Autoregression. #' #' @examples -#' \dontrun{ -#' data(us_fiscal_lsuw) +#' data(ilo_cubic_panel) #' spec = specify_bvarPANEL$new( -#' data = us_fiscal_lsuw, +#' data = ilo_cubic_panel, #' p = 4 #' ) -#' } #' #' @export specify_bvarPANEL = R6::R6Class( @@ -402,14 +409,13 @@ specify_bvarPANEL = R6::R6Class( #' Returns the data matrices as the DataMatricesBVARPANEL object. #' #' @examples - #' \dontrun{ - #' data(us_fiscal_lsuw) - #' spec = specify_bsvar$new( - #' data = us_fiscal_lsuw, + #' data(ilo_cubic_panel) + #' spec = specify_bvarPANEL$new( + #' data = ilo_cubic_panel, #' p = 4 #' ) #' spec$get_data_matrices() - #' } + #' get_data_matrices = function() { self$data_matrices$clone() }, # END get_data_matrices @@ -418,14 +424,13 @@ specify_bvarPANEL = R6::R6Class( #' Returns the prior specification as the PriorBVARPANEL object. #' #' @examples - #' \dontrun{ - #' data(us_fiscal_lsuw) - #' spec = specify_bsvar$new( - #' data = us_fiscal_lsuw, + #' data(ilo_cubic_panel) + #' spec = specify_bvarPANEL$new( + #' data = ilo_cubic_panel, #' p = 4 #' ) #' spec$get_prior() - #' } + #' get_prior = function() { self$prior$clone() }, # END get_prior @@ -434,14 +439,13 @@ specify_bvarPANEL = R6::R6Class( #' Returns the starting values as the StartingValuesBVARPANEL object. #' #' @examples - #' \dontrun{ - #' data(us_fiscal_lsuw) - #' spec = specify_bsvar$new( - #' data = us_fiscal_lsuw, + #' data(ilo_cubic_panel) + #' spec = specify_bvarPANEL$new( + #' data = ilo_cubic_panel, #' p = 4 #' ) #' spec$get_starting_values() - #' } + #' get_starting_values = function() { self$starting_values$clone() } # END get_starting_values @@ -463,14 +467,17 @@ specify_bvarPANEL = R6::R6Class( #' @seealso \code{\link{specify_bvarPANEL}} #' #' @examples -#' \dontrun{ #' # This is a function that is used within estimate() -#' data(us_fiscal_lsuw) -#' specification = specify_bsvar$new(us_fiscal_lsuw, p = 4) +#' data(ilo_cubic_panel) #' set.seed(123) -#' estimate = estimate(specification, 50) -#' class(estimate) -#' } +#' specification = specify_bvarPANEL$new( +#' data = ilo_cubic_panel, +#' p = 4 +#' ) +#' +#' posterior = estimate(specification, 50) +#' class(posterior) +#' #' @export specify_posterior_bvarPANEL = R6::R6Class( "PosteriorBVARPANEL", @@ -508,13 +515,16 @@ specify_posterior_bvarPANEL = R6::R6Class( #' Returns a list containing Bayesian estimation output. #' #' @examples - #' \dontrun{ - #' data(us_fiscal_lsuw) - #' specification = specify_bsvar$new(us_fiscal_lsuw) + #' data(ilo_cubic_panel) #' set.seed(123) - #' estimate = estimate(specification, 50) - #' estimate$get_posterior() - #' } + #' specification = specify_bvarPANEL$new( + #' data = ilo_cubic_panel, + #' p = 4 + #' ) + #' + #' posterior = estimate(specification, 50) + #' posterior$get_posterior() + #' get_posterior = function(){ self$posterior }, # END get_posterior @@ -525,19 +535,19 @@ specify_posterior_bvarPANEL = R6::R6Class( #' MCMC estimation using \code{estimate()}. #' #' @examples - #' \dontrun{ - #' data(us_fiscal_lsuw) - #' - #' # specify the model and set seed - #' specification = specify_bsvar$new(us_fiscal_lsuw, p = 4) + #' data(ilo_cubic_panel) #' set.seed(123) + #' specification = specify_bvarPANEL$new( + #' data = ilo_cubic_panel, + #' p = 4 + #' ) #' #' # run the burn-in #' burn_in = estimate(specification, 10) #' #' # estimate the model #' posterior = estimate(burn_in, 10) - #' } + #' get_last_draw = function(){ self$last_draw$clone() } # END get_last_draw diff --git a/man/specify_bvarPANEL.Rd b/man/specify_bvarPANEL.Rd index 8d7e888..99f1ff0 100644 --- a/man/specify_bvarPANEL.Rd +++ b/man/specify_bvarPANEL.Rd @@ -8,53 +8,48 @@ The class BVARPANEL presents complete specification for the Bayesian Panel Vector Autoregression. } \examples{ -\dontrun{ -data(us_fiscal_lsuw) +data(ilo_cubic_panel) spec = specify_bvarPANEL$new( - data = us_fiscal_lsuw, + data = ilo_cubic_panel, p = 4 ) -} ## ------------------------------------------------ ## Method `specify_bvarPANEL$get_data_matrices` ## ------------------------------------------------ -\dontrun{ -data(us_fiscal_lsuw) -spec = specify_bsvar$new( - data = us_fiscal_lsuw, +data(ilo_cubic_panel) +spec = specify_bvarPANEL$new( + data = ilo_cubic_panel, p = 4 ) spec$get_data_matrices() -} + ## ------------------------------------------------ ## Method `specify_bvarPANEL$get_prior` ## ------------------------------------------------ -\dontrun{ -data(us_fiscal_lsuw) -spec = specify_bsvar$new( - data = us_fiscal_lsuw, +data(ilo_cubic_panel) +spec = specify_bvarPANEL$new( + data = ilo_cubic_panel, p = 4 ) spec$get_prior() -} + ## ------------------------------------------------ ## Method `specify_bvarPANEL$get_starting_values` ## ------------------------------------------------ -\dontrun{ -data(us_fiscal_lsuw) -spec = specify_bsvar$new( - data = us_fiscal_lsuw, +data(ilo_cubic_panel) +spec = specify_bvarPANEL$new( + data = ilo_cubic_panel, p = 4 ) spec$get_starting_values() -} + } \section{Public fields}{ \if{html}{\out{
}} @@ -113,14 +108,13 @@ Returns the data matrices as the DataMatricesBVARPANEL object. \subsection{Examples}{ \if{html}{\out{
}} -\preformatted{\dontrun{ -data(us_fiscal_lsuw) -spec = specify_bsvar$new( - data = us_fiscal_lsuw, +\preformatted{data(ilo_cubic_panel) +spec = specify_bvarPANEL$new( + data = ilo_cubic_panel, p = 4 ) spec$get_data_matrices() -} + } \if{html}{\out{
}} @@ -138,14 +132,13 @@ Returns the prior specification as the PriorBVARPANEL object. \subsection{Examples}{ \if{html}{\out{
}} -\preformatted{\dontrun{ -data(us_fiscal_lsuw) -spec = specify_bsvar$new( - data = us_fiscal_lsuw, +\preformatted{data(ilo_cubic_panel) +spec = specify_bvarPANEL$new( + data = ilo_cubic_panel, p = 4 ) spec$get_prior() -} + } \if{html}{\out{
}} @@ -163,14 +156,13 @@ Returns the starting values as the StartingValuesBVARPANEL object. \subsection{Examples}{ \if{html}{\out{
}} -\preformatted{\dontrun{ -data(us_fiscal_lsuw) -spec = specify_bsvar$new( - data = us_fiscal_lsuw, +\preformatted{data(ilo_cubic_panel) +spec = specify_bvarPANEL$new( + data = ilo_cubic_panel, p = 4 ) spec$get_starting_values() -} + } \if{html}{\out{
}} diff --git a/man/specify_panel_data_matrices.Rd b/man/specify_panel_data_matrices.Rd index 284a246..219f061 100644 --- a/man/specify_panel_data_matrices.Rd +++ b/man/specify_panel_data_matrices.Rd @@ -7,6 +7,21 @@ The class DataMatricesBVARPANEL presents the data matrices of dependent variables, \eqn{\mathbf{Y}_c}, and regressors, \eqn{\mathbf{X}_c}, for the Bayesian Panel VAR model for all countries \eqn{c = 1, ..., C}. +} +\examples{ +data(ilo_cubic_panel) +YX = specify_panel_data_matrices$new(data = ilo_cubic_panel, p = 4) +length(YX$Y); names(YX$Y) + + +## ------------------------------------------------ +## Method `specify_panel_data_matrices$get_data_matrices` +## ------------------------------------------------ + +data(ilo_cubic_panel) +YX = specify_panel_data_matrices$new(ilo_cubic_panel) +YX$get_data_matrices() + } \section{Public fields}{ \if{html}{\out{
}} @@ -59,6 +74,17 @@ Returns the data matrices DataMatricesBVARPANEL as a \code{list}. \if{html}{\out{
}}\preformatted{specify_panel_data_matrices$get_data_matrices()}\if{html}{\out{
}} } +\subsection{Examples}{ +\if{html}{\out{
}} +\preformatted{data(ilo_cubic_panel) +YX = specify_panel_data_matrices$new(ilo_cubic_panel) +YX$get_data_matrices() + +} +\if{html}{\out{
}} + +} + } \if{html}{\out{
}} \if{html}{\out{}} diff --git a/man/specify_posterior_bvarPANEL.Rd b/man/specify_posterior_bvarPANEL.Rd index 4a2ebe2..4fa0d84 100644 --- a/man/specify_posterior_bvarPANEL.Rd +++ b/man/specify_posterior_bvarPANEL.Rd @@ -11,44 +11,50 @@ Note that due to the thinning of the MCMC output the starting value in element element \code{posterior}. } \examples{ -\dontrun{ # This is a function that is used within estimate() -data(us_fiscal_lsuw) -specification = specify_bsvar$new(us_fiscal_lsuw, p = 4) +data(ilo_cubic_panel) set.seed(123) -estimate = estimate(specification, 50) -class(estimate) -} +specification = specify_bvarPANEL$new( + data = ilo_cubic_panel, + p = 4 +) + +posterior = estimate(specification, 50) +class(posterior) + ## ------------------------------------------------ ## Method `specify_posterior_bvarPANEL$get_posterior` ## ------------------------------------------------ -\dontrun{ -data(us_fiscal_lsuw) -specification = specify_bsvar$new(us_fiscal_lsuw) +data(ilo_cubic_panel) set.seed(123) -estimate = estimate(specification, 50) -estimate$get_posterior() -} +specification = specify_bvarPANEL$new( + data = ilo_cubic_panel, + p = 4 +) + +posterior = estimate(specification, 50) +posterior$get_posterior() + ## ------------------------------------------------ ## Method `specify_posterior_bvarPANEL$get_last_draw` ## ------------------------------------------------ -\dontrun{ -data(us_fiscal_lsuw) - -# specify the model and set seed -specification = specify_bsvar$new(us_fiscal_lsuw, p = 4) +data(ilo_cubic_panel) set.seed(123) +specification = specify_bvarPANEL$new( + data = ilo_cubic_panel, + p = 4 +) # run the burn-in burn_in = estimate(specification, 10) # estimate the model posterior = estimate(burn_in, 10) -} + } \seealso{ \code{\link{specify_bvarPANEL}} @@ -107,13 +113,16 @@ Returns a list containing Bayesian estimation output. \subsection{Examples}{ \if{html}{\out{
}} -\preformatted{\dontrun{ -data(us_fiscal_lsuw) -specification = specify_bsvar$new(us_fiscal_lsuw) +\preformatted{data(ilo_cubic_panel) set.seed(123) -estimate = estimate(specification, 50) -estimate$get_posterior() -} +specification = specify_bvarPANEL$new( + data = ilo_cubic_panel, + p = 4 +) + +posterior = estimate(specification, 50) +posterior$get_posterior() + } \if{html}{\out{
}} @@ -133,19 +142,19 @@ MCMC estimation using \code{estimate()}. \subsection{Examples}{ \if{html}{\out{
}} -\preformatted{\dontrun{ -data(us_fiscal_lsuw) - -# specify the model and set seed -specification = specify_bsvar$new(us_fiscal_lsuw, p = 4) +\preformatted{data(ilo_cubic_panel) set.seed(123) +specification = specify_bvarPANEL$new( + data = ilo_cubic_panel, + p = 4 +) # run the burn-in burn_in = estimate(specification, 10) # estimate the model posterior = estimate(burn_in, 10) -} + } \if{html}{\out{
}} From 4777241f5e5b470125ebc4a1c3bb3a4a0a64800e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sat, 11 May 2024 23:31:04 +1000 Subject: [PATCH 30/32] full examples included #3 #2 --- R/bvarPANELs-package.R | 6 ++++++ R/estimate.bvarPANEL.R | 7 +++++++ man/bvarPANELs-package.Rd | 21 +++++++++++++++++++++ man/estimate.BVARPANEL.Rd | 8 ++++++++ man/estimate.PosteriorBVARPANEL.Rd | 8 ++++++++ 5 files changed, 50 insertions(+) diff --git a/R/bvarPANELs-package.R b/R/bvarPANELs-package.R index 5dd105a..3131f1a 100644 --- a/R/bvarPANELs-package.R +++ b/R/bvarPANELs-package.R @@ -38,4 +38,10 @@ #' @note This package is currently in active development. #' @author Tomasz Woźniak \email{wozniak.tom@pm.me} #' @keywords package models ts +#' #' @examples +#' data(ilo_cubic_panel) # load the data +#' set.seed(123) +#' specification = specify_bvarPANEL$new(ilo_cubic_panel) # specify the model +#' burn_in = estimate(specification, 10) # run the burn-in +#' posterior = estimate(burn_in, 10) # estimate the model "_PACKAGE" diff --git a/R/estimate.bvarPANEL.R b/R/estimate.bvarPANEL.R index 6f58991..ddb0b43 100644 --- a/R/estimate.bvarPANEL.R +++ b/R/estimate.bvarPANEL.R @@ -36,6 +36,13 @@ #' #' @method estimate BVARPANEL #' +#' @examples +#' data(ilo_cubic_panel) # load the data +#' set.seed(123) +#' specification = specify_bvarPANEL$new(ilo_cubic_panel) # specify the model +#' burn_in = estimate(specification, 10) # run the burn-in +#' posterior = estimate(burn_in, 10) # estimate the model +#' #' @export estimate.BVARPANEL <- function( specification, diff --git a/man/bvarPANELs-package.Rd b/man/bvarPANELs-package.Rd index dfa56f2..6df853d 100644 --- a/man/bvarPANELs-package.Rd +++ b/man/bvarPANELs-package.Rd @@ -16,6 +16,27 @@ This package is currently in active development. \author{ Tomasz Woźniak \email{wozniak.tom@pm.me} } +\keyword{#} +\keyword{#'} +\keyword{10)} +\keyword{=} +\keyword{@examples} +\keyword{burn-in} +\keyword{burn_in} +\keyword{data} +\keyword{data(ilo_cubic_panel)} +\keyword{estimate} +\keyword{estimate(burn_in,} +\keyword{estimate(specification,} +\keyword{load} +\keyword{model} \keyword{models} \keyword{package} +\keyword{posterior} +\keyword{run} +\keyword{set.seed(123)} +\keyword{specification} +\keyword{specify} +\keyword{specify_bvarPANEL$new(ilo_cubic_panel)} +\keyword{the} \keyword{ts} diff --git a/man/estimate.BVARPANEL.Rd b/man/estimate.BVARPANEL.Rd index 111824f..8d063c4 100644 --- a/man/estimate.BVARPANEL.Rd +++ b/man/estimate.BVARPANEL.Rd @@ -40,6 +40,14 @@ country \code{c}, \eqn{X_c} is a \code{T_c x K} matrix of explanatory variables, \eqn{E_c} is an \code{T_c x N} matrix of error terms, and \eqn{A_c} is an \code{NxK} matrix of country-specific autoregressive slope coefficients and parameters on deterministic terms in \eqn{X_c}. +} +\examples{ +data(ilo_cubic_panel) # load the data +set.seed(123) +specification = specify_bvarPANEL$new(ilo_cubic_panel) # specify the model +burn_in = estimate(specification, 10) # run the burn-in +posterior = estimate(burn_in, 10) # estimate the model + } \seealso{ \code{\link{specify_bvarPANEL}}, \code{\link{specify_posterior_bvarPANEL}} diff --git a/man/estimate.PosteriorBVARPANEL.Rd b/man/estimate.PosteriorBVARPANEL.Rd index 7bad6df..2b65bb3 100644 --- a/man/estimate.PosteriorBVARPANEL.Rd +++ b/man/estimate.PosteriorBVARPANEL.Rd @@ -42,6 +42,14 @@ country \code{c}, \eqn{X_c} is a \code{T_c x K} matrix of explanatory variables, \eqn{E_c} is an \code{T_c x N} matrix of error terms, and \eqn{A_c} is an \code{NxK} matrix of country-specific autoregressive slope coefficients and parameters on deterministic terms in \eqn{X_c}. +} +\examples{ +data(ilo_cubic_panel) # load the data +set.seed(123) +specification = specify_bvarPANEL$new(ilo_cubic_panel) # specify the model +burn_in = estimate(specification, 10) # run the burn-in +posterior = estimate(burn_in, 10) # estimate the model + } \seealso{ \code{\link{specify_bvarPANEL}}, \code{\link{specify_posterior_bvarPANEL}} From 692f293f5df387d097fab8b57ba2e87ef8b02038 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sat, 11 May 2024 23:37:23 +1000 Subject: [PATCH 31/32] Create NEWS.md --- NEWS.md | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 NEWS.md diff --git a/NEWS.md b/NEWS.md new file mode 100644 index 0000000..c5bbd76 --- /dev/null +++ b/NEWS.md @@ -0,0 +1,5 @@ +# bvarPANELs 0.0.1.9000 + +1. This package is alive! [#1](https://github.com/bsvars/bvarPANELs/issues/1) +2. The package includes the Gibbs sampler for a Bayesian Hierarchical Panel Vector Autoregression model. [#2](https://github.com/bsvars/bvarPANELs/issues/2) +3. The package includes labour market dynamic panel data. [#3](https://github.com/bsvars/bvarPANELs/issues/3) From b7a7d34ffe30b41030353dd3d7033c4dc298b6ba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Sat, 11 May 2024 23:58:08 +1000 Subject: [PATCH 32/32] included tests for Gibbs sampler #2 --- DESCRIPTION | 6 +++- inst/tinytest/test_estimate_bvarpanel.R | 44 +++++++++++++++++++++++++ tests/tinytest.R | 7 ++++ 3 files changed, 56 insertions(+), 1 deletion(-) create mode 100644 inst/tinytest/test_estimate_bvarpanel.R create mode 100644 tests/tinytest.R diff --git a/DESCRIPTION b/DESCRIPTION index 171bf5e..e364d28 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -4,7 +4,10 @@ Title: Forecasting with Bayesian Hierarchical Panel Vector Autoregressions Description: Forecasting a multi-country time series panel data using Bayesian Vector Autoregressions with a three-level country-global hierarchical prior structure. Version: 0.0.1.9000 Date: 2024-04-27 -Author: c(person(given="Tomasz", family="Woźniak", email="wozniak.tom@pm.me", role = c("aut", "cre"), comment = c(ORCID = "0000-0003-2212-2378")), person(given="Miguel", family="Sanchez-Martinez", role = "ctb"), person(given="International Labour Organization", role = "cph")) +Author: c(person(given="Tomasz", family="Woźniak", email="wozniak.tom@pm.me", role = + c("aut", "cre"), comment = c(ORCID = "0000-0003-2212-2378")), + person(given="Miguel", family="Sanchez-Martinez", role = "ctb"), + person(given="International Labour Organization", role = "cph")) Maintainer: Tomasz Woźniak License: GPL (>= 3) Depends: @@ -16,6 +19,7 @@ Imports: RcppProgress, RcppTN LinkingTo: bsvars, Rcpp, RcppArmadillo, RcppProgress, RcppTN +Suggests: tinytest LazyData: true Encoding: UTF-8 RoxygenNote: 7.3.1 diff --git a/inst/tinytest/test_estimate_bvarpanel.R b/inst/tinytest/test_estimate_bvarpanel.R new file mode 100644 index 0000000..3dd8fdd --- /dev/null +++ b/inst/tinytest/test_estimate_bvarpanel.R @@ -0,0 +1,44 @@ +data(ilo_cubic_panel) + +set.seed(1) +specification_no1 <- specify_bvarPANEL$new(ilo_cubic_panel) +run_no1 <- estimate(specification_no1, 3, 1, show_progress = FALSE) + +set.seed(1) +specification_no2 <- specify_bvarPANEL$new(ilo_cubic_panel) +run_no2 <- estimate(specification_no2, 3, 1, show_progress = FALSE) + +set.seed(1) +run_no3 <- ilo_cubic_panel |> + specify_bvarPANEL$new() |> + estimate(S = 3, thin = 1, show_progress = FALSE) + +expect_identical( + run_no1$last_draw$starting_values$A[1,1], + run_no2$last_draw$starting_values$A[1,1], + info = "estimate_bvarPANEL: the last_draw(s) of two runs to be identical." +) + +expect_identical( + run_no1$posterior$A[1,1,1], + run_no2$posterior$A[1,1,1], + info = "estimate_bvarPANEL: the first draws of two runs to be identical." +) + +expect_identical( + run_no1$last_draw$starting_values$A[1,1], + run_no3$last_draw$starting_values$A[1,1], + info = "estimate_bvarPANEL: the last_draw(s) of a normal and pipe run to be identical." +) + +# a test of a good setting of S and thin +expect_error( + estimate(specification_no1, 3, 2, show_progress = FALSE), + info = "Argument S is not a positive integer multiplication of argument thin." +) + +expect_error( + estimate(specification_no1, 2, 3, show_progress = FALSE), + info = "Argument S is not a positive integer multiplication of argument thin." +) + diff --git a/tests/tinytest.R b/tests/tinytest.R new file mode 100644 index 0000000..b1418c8 --- /dev/null +++ b/tests/tinytest.R @@ -0,0 +1,7 @@ + +if ( requireNamespace("tinytest", quietly=TRUE) ){ + # To skip running the tests on CRAN make them run only on the developer version + home <- length(unclass(packageVersion("bsvars"))[[1]]) == 4 + tinytest::test_package("bvarPANELs") +} +