Данные
Здесь я подготовил и визуализировал данные.
# Load packages
library(tidyverse)
library(rjags)
# Set seed for reproduciblility
set.seed(199)
mkt.bull <- rnorm(150, 2, 1.5)
mkt.bear <- rnorm(150, -1, 2.5)
x <- c(mkt.bear,mkt.bull)
portfolio.bull <- rnorm(150, 1.75, 1.6)
portfolio.bear <- rnorm(150, -0.5, 2.3)
y <- c(portfolio.bear,portfolio.bull)
# Create example data frame
dat <- data.frame(x = x, y = y, regime = c(rep("bear", 150), rep("bull", 150)),
stringsAsFactors = FALSE)
# Plot the sample distribution
dat$regime <- factor(dat$regime, levels = c("bear", "bull"))
# Create a plot
ggplot(dat, aes(x = y, color = regime)) +
geom_density()
Есть два режима: bear и bull. y для этих режимов нормально распределены. Похоже, что ОП хочет оценить среднее и стандартное отклонение y, обусловленное этими состояниями.
Максимальное правдоподобие
Вот один из способов использования максимальной вероятности для оценки параметров с использованием пакета stats4.
# Load the infer package
library(stats4)
# Split the data
y_bull <- dat %>% filter(regime %in% "bull") %>% pull("y")
y_bear <- dat %>% filter(regime %in% "bear") %>% pull("y")
# Define the log-likelihood function
LogLike_bull <- function(Mean, Sigma){
R <- suppressWarnings(dnorm(y_bull, Mean, Sigma))
return(-sum(log(R)))
}
LogLike_bear <- function(Mean, Sigma){
R <- suppressWarnings(dnorm(y_bear, Mean, Sigma))
return(-sum(log(R)))
}
mle(minuslogl = LogLike_bull, start = list(Mean = 1, Sigma = 1))
# Call:
# mle(minuslogl = LogLike_bull, start = list(Mean = 1, Sigma = 1))
#
# Coefficients:
# Mean Sigma
# 1.703099 1.482619
mle(minuslogl = LogLike_bear, start = list(Mean = 1, Sigma = 1))
# Call:
# mle(minuslogl = LogLike_bear, start = list(Mean = 1, Sigma = 1))
#
# Coefficients:
# Mean Sigma
# -0.616106 2.340852
Параметры для bull являются средними = 1.703 и стандартным отклонением = 1.483. Параметры для bear - это среднее значение = -0.616 и стандартное отклонение = 2.341. Они близки к истинным ценностям.
Байсеевский анализ
Вот попытка использовать байесовский анализ для решения этого вопроса с помощью jags и пакета rjags.
Я использовал байесовскую модель для оценки alpha (среднее значение y на bear), beta (разница y на bear и bull) и сигмы (стандарт отклонение y на bear и bull) с использованием 10000 итераций.
# Define the Bayesian model
model <- "model{
for(i in 1:length(Y)) {
Y[i] ~ dnorm(Mean[i], s[X[i]]^(-2))
Mean[i] <- alpha + beta[X[i]]
}
alpha ~ dnorm(0, 5^(-2))
beta[1] <- 0
beta[2] ~ dnorm(0, 5^(-2))
s[1] ~ dunif(0, 10)
s[2] ~ dunif(0, 10)
}"
# Compile the model
jags_model <- jags.model(
textConnection(model),
data = list(Y = dat$y, X = dat$regime),
n.chains = 3,
inits = list(.RNG.name = "base::Wichmann-Hill", .RNG.seed = 10)
)
# Simulate the posterior
jags_sim <- coda.samples(model = jags_model,
variable.names = c("alpha", "beta", "s"),
n.iter = 10000)
# Plot the posterior
plot(jags_sim)
График показывает, что оценки хорошо перемешаны.
# See the summary
summary(jags_sim)
Iterations = 1001:11000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 10000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
alpha -0.614 0.19436 0.0011222 0.0027201
beta[1] 0.000 0.00000 0.0000000 0.0000000
beta[2] 2.315 0.23099 0.0013336 0.0032666
s[1] 2.369 0.13768 0.0007949 0.0010393
s[2] 1.500 0.08836 0.0005102 0.0006727
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
alpha -0.9838 -0.7471 -0.6147 -0.4857 -0.2225
beta[1] 0.0000 0.0000 0.0000 0.0000 0.0000
beta[2] 1.8582 2.1622 2.3174 2.4760 2.7586
s[1] 2.1225 2.2722 2.3632 2.4564 2.6611
s[2] 1.3368 1.4390 1.4959 1.5579 1.6813
Среднее значение alpha для bear равно -0.614, что аналогично фактическому значению -0.5. Среднее значение beta[2] равно 2.315. Если мы добавим alpha и beta[2], мы получим 1.701, тогда как фактическое значение будет 1.7. Мы также получили s[1] и s[2] как 2.369 и 1.5, которые аналогичны 2.3 и 1.6 соответственно.
Бутстрапирование
Вот еще один подход к использованию начальной загрузки для оценки alpha, beta и стандартного отклонения, основанный на пакете infer.
# Load the infer package
library(infer)
set.seed(199)
# Split the data
dat_bull <- dat %>% filter(regime %in% "bull")
dat_bear <- dat %>% filter(regime %in% "bear")
# Calcualte the values in bull
dat_bull2 <- dat_bull %>%
# Specify the response variable
specify(response = y) %>%
# Generate 10000 bootstrap samples
generate(reps = 10000, type = "bootstrap")
summary_bull <- dat_bull2 %>%
summarise(mean_y = mean(y), sd_y = sd(y))
# Calcualte the values in bear
dat_ear2 <- dat_bear %>%
# Specify the response variable
specify(response = y) %>%
# Generate 10000 bootstrap samples
generate(reps = 10000, type = "bootstrap")
summary_bear <- dat_ear2 %>%
summarise(mean_y = mean(y), sd_y = sd(y))
Теперь мы можем напечатать результаты. Все они похожи на истинные значения.
# The mean of bull
mean(summary_bull$mean_y)
# [1] 1.702693
# The standard deviation of bear
mean(summary_bull$sd_y)
# [1] 1.480158
# The mean of bear
mean(summary_bear$mean_y)
# [1] -0.6165585
# The standard deviation of bear
mean(summary_bear$sd_y)
# [1] 2.337042