时间戳

首页 > 解决方案 > 使用回归参数作为 rnorm 的平均值

r - 使用回归参数作为 rnorm 的平均值

问题描述

我想测试一个模型，其中假设正常的随机变量的分布取决于另一个随机变量的状态，该状态根据马尔可夫链切换状态。第一步是：

假设简单的线性模型：

lm(y~x, data=data)

假设 x 切换状态，我想找到分布的参数。

例如：

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)

在上面的示例中，x 可以建模为msmFit具有两种状态的马尔可夫切换模型 ( )，一种牛市和一种熊市。而不是用lm来解决问题，

lm(y~x)

由于这两个系列显然是非线性的，我想运行一个回归，其中参数取决于制度。这可以通过最大似然来完成，但第一步是将 y 的分布定义为：

y_i | x, S_t ~ N(alpha + beta_{i,s_t}); sigma^2)

如何编码上述公式？我想这不能使用 rnorm 来完成。还有其他方法吗？

谢谢

标签： rregression

解决方案

数据

在这里，我准备并可视化了数据。

# 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这些制度，两者都是正态分布的。OP 似乎想要估计以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包。

我运行了一个贝叶斯模型，使用 10000 次迭代来估计（ on的alpha平均值）、（on和的差）和 sigma（on和的标准差）。ybearbetaybearbullybearbull

# 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

alphaon的平均值为bear，-0.614与实际值相近-0.5。的平均值beta[2]是2.315。如果我们加上alphaand beta[2]，我们得到了1.701，而实际值为1.7。我们还得到了s[1]和s[2]作为2.369和1.5，它们分别类似于2.3和1.6。

自举

这是使用 bootstrap 估计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

推荐阅读