I would like if this code is coherent and provide correct plots?

Load necessary libraries

library(MCMCpack)
library(MASS)
library(ggplot2)
library(dplyr)
library(tidyr)
library(imputeTS)

Load the user data

data <- read.csv("C:/Users/TBS/Documents/Data-Was.csv")

Define output directory

output_dir <- "C:/Users/TBS/Documents/Output/"
if (!dir.exists(output_dir)) {
dir.create(output_dir, recursive = TRUE)
}

Prepare country-specific data

countries <- unique(data$Country)
N <- length(countries) # Number of units (countries)
time_column <- "Date"
T <- length(unique(data[[time_column]])) # Number of time points
M <- 3 # Variables: TempAnom, CPI, GDP
K <- 2 # Two regimes
n_iter <- 5000 # Number of iterations for Gibbs sampler
burnin <- 1000 # Burn-in period for Gibbs sampling

Initialize matrices for TempAnom, CPI, and GDP

TempAnom <- matrix(NA, nrow = N, ncol = T)
CPI <- matrix(NA, nrow = N, ncol = T)
GDP <- matrix(NA, nrow = N, ncol = T)

Fill matrices with data for each country

for (i in 1:N) {
country_data <- data %>%
filter(Country == countries[i]) %>%
arrange(!!sym(time_column))

if (nrow(country_data) != T) {
stop(paste("Data for country", countries[i], "has inconsistent time periods."))
}

Replace non-positive values with a small positive value

country_data$TempAnom <- ifelse(country_data$TempAnom <= 0, 0.01, country_data$TempAnom)
country_data$CPI <- ifelse(country_data$CPI <= 0, 0.01, country_data$CPI)
country_data$GDP <- ifelse(country_data$GDP <= 0, 0.01, country_data$GDP)

Logarithmic transformation

TempAnom[i, ] <- log(na_interpolation(country_data$TempAnom))
CPI[i, ] <- log(na_interpolation(country_data$CPI))
GDP[i, ] <- log(na_interpolation(country_data$GDP))
}

Create a 3-dimensional array for data

Y <- array(NA, dim = c(N, T, M))
Y[, , 1] <- TempAnom
Y[, , 2] <- CPI
Y[, , 3] <- GDP

FFBS function

ffbs <- function(log_probs, P, T) {
K <- ncol(P)
S <- numeric(T)
alpha <- matrix(0, nrow = T, ncol = K)

alpha[1, ] <- log_probs[1, ] - log(sum(exp(log_probs[1, ])))
for (t in 2:T) {
for (k in 1:K) {
alpha[t, k] <- log_probs[t, k] + log(sum(exp(alpha[t-1, ] + log(P[, k]))))
}
}

S[T] <- sample(1:K, 1, prob = exp(alpha[T, ]))
for (t in (T-1):1) {
trans_probs <- exp(alpha[t, ] + log(P[, S[t+1]]))
S[t] <- sample(1:K, 1, prob = trans_probs / sum(trans_probs))
}
return(S)
}

Gibbs sampler function (updated with corrections)

gibbs_sampler_block <- function(Y, N, T, M, K, iterations, burnin) {
beta <- array(0, dim = c(iterations, N, K, M + 1)) # Include intercept term
Sigma <- array(0, dim = c(iterations, N, K, M + 1, M + 1)) # Covariance matrices
alpha <- array(0, dim = c(N, K, K, 2)) # Transition parameters
P <- array(0, dim = c(N, K, K)) # Transition probability matrix

for (iter in 2:iterations) {
for (i in 1:N) { # Loop over countries
country_data <- Y[i, , ]

  for (k in 1:K) {  # Loop over regimes
    # Design matrix with intercept
    X_ik <- cbind(1, country_data[, 1:M])
    y_ik <- country_data[, k]
    
    # Posterior for beta
    XtX <- t(X_ik) %*% X_ik + diag(1e-3, M + 1)  # Regularization
    XtY <- t(X_ik) %*% y_ik
    beta_mean <- solve(XtX) %*% XtY
    beta_var <- solve(XtX)
    beta[iter, i, k, ] <- mvrnorm(1, beta_mean, beta_var)
    
    # Residuals as a matrix
    residuals <- y_ik - X_ik %*% beta[iter, i, k, ]
    residuals <- matrix(residuals, nrow = T, ncol = 1)  # Ensure residuals are a T x 1 matrix
    
    # Update covariance matrix (Sigma)
    S <- (t(X_ik) %*% diag(residuals[, 1]^2) %*% X_ik) / T + diag(0.001, M + 1)
    Sigma[iter, i, k, , ] <- riwish(T + M, S)
  }
  
  # Transition probabilities with centered parametrization
  for (l in 1:K) {  # Previous state
    denom <- 0
    for (k in 1:K) {  # Current state
      if (k != K) {  # Reference state
        Vt <- apply(country_data, 2, mean)  # Average over time
        ci <- mean(Vt)  # Threshold parameter
        log_prob <- sum((Vt - ci) * alpha[i, k, l, ])
        P[i, k, l] <- exp(log_prob)  # Numerator
        denom <- denom + exp(log_prob)
      } else {
        P[i, k, l] <- 1  # Reference state
        denom <- denom + 1
      }
    }
    # Normalize probabilities
    for (k in 1:K) {
      P[i, k, l] <- P[i, k, l] / denom
    }
  }
}

}

return(list(beta = beta, Sigma = Sigma, P = P, alpha = alpha))
}

Run the Gibbs sampler

posterior_samples <- gibbs_sampler_block(Y, N, T, M, K, n_iter, burnin)

Extract posterior samples

beta_samples <- posterior_samples$beta
Sigma_samples <- posterior_samples$Sigma
P_samples <- posterior_samples$P
alpha_samples <- posterior_samples$alpha

Plotting the results

beta_df <- data.frame()
sigma_df <- data.frame()

for (i in 1:N) {
country_name <- countries[i]
for (k in 1:K) {
for (iter in (burnin + 1):n_iter) {
beta_df <- rbind(beta_df, data.frame(
country = country_name,
Regime = paste("Regime", k),
Variable = "TempAnom-CPI",
Intercept = beta_samples[iter, i, k, 2]
))

  beta_df <- rbind(beta_df, data.frame(
    country = country_name,
    Regime = paste("Regime", k),
    Variable = "TempAnom-GDP",
    Intercept = beta_samples[iter, i, k, 3]
  ))
  
  sigma_df <- rbind(sigma_df, data.frame(
    country = country_name,
    Regime = paste("Regime", k),
    Variable = "TempAnom-CPI",
    Volatility = sqrt(Sigma_samples[iter, i, k, 2, 2])
  ))
  
  sigma_df <- rbind(sigma_df, data.frame(
    country = country_name,
    Regime = paste("Regime", k),
    Variable = "TempAnom-GDP",
    Volatility = sqrt(Sigma_samples[iter, i, k, 3, 3])
  ))
}

}
}

Save plots for each country

for (country in countries) {
beta_data <- beta_df %>% filter(country == country)
sigma_data <- sigma_df %>% filter(country == country)

interceptplot <- ggplot(beta_data, aes(x = Intercept, color = Regime)) +
geom_density(adjust = 1.5, alpha = 0.5) +
facet_wrap(~Variable, scales = "free") +
labs(title = paste("Intercepts for", country), x = "Intercept", y = "Density") +
theme_minimal()
ggsave(paste0(output_dir, "Intercepts", country, ".png"), intercept_plot, width = 10, height = 5)

volatilityplot <- ggplot(sigma_data, aes(x = Volatility, color = Regime)) +
geom_density(adjust = 1.5, alpha = 0.5) +
facet_wrap(~Variable, scales = "free") +
labs(title = paste("Volatilities for", country), x = "Volatility", y = "Density") +
theme_minimal()
ggsave(paste0(output_dir, "Volatilities", country, ".png"), volatility_plot, width = 10, height = 5)
}