11  时间序列回归

library(cmdstanr)
library(zoo)
library(xts) # xts 依赖 zoo
library(fGarch)
library(INLA)
library(mgcv)
library(tensorflow)
library(ggplot2)
library(bayesplot)

11.1 随机波动率模型

# 美团上市至 2023-07-15
library(zoo)
library(xts)
library(ggplot2)
theme_classic() +
labs(x = "日期", y = "股价")

$\text{对数收益率} = \ln(\text{今日收盘价} / \text{昨日收盘价} ) = \ln (1 + \text{普通收益率})$

meituan_log_return <- diff(log(meituan[, "3690.HK.Adjusted"]))[-1]
autoplot(meituan_log_return) +
theme_classic() +
labs(x = "日期", y = "对数收益率")
ggplot(data = meituan_log_return, aes(x = 3690.HK.Adjusted)) +
geom_histogram(color = "black", fill = "gray", bins = 30) +
theme_classic() +
labs(x = "对数收益率", y = "频数（天数）")

acf(meituan_log_return, main = "")

Box.test(meituan_log_return, lag = 12, type = "Ljung")
#>
#>  Box-Ljung test
#>
#> data:  meituan_log_return
#> X-squared = 35.669, df = 12, p-value = 0.0003661

# ARCH 效应的检验
Box.test((meituan_log_return - mean(meituan_log_return))^2,
lag = 12, type = "Ljung")
#>
#>  Box-Ljung test
#>
#> data:  (meituan_log_return - mean(meituan_log_return))^2
#> X-squared = 124.54, df = 12, p-value < 2.2e-16

11.1.1 Stan 框架

\begin{aligned} y_t &= \epsilon_t \exp(h_t / 2) \\ h_{t+1} &= \mu + \phi (h_t - \mu) + \delta_t \sigma \\ h_1 &\sim \textsf{normal}\left( \mu, \frac{\sigma}{\sqrt{1 - \phi^2}} \right) \\ \epsilon_t &\sim \textsf{normal}(0,1) \\ \delta_t &\sim \textsf{normal}(0,1) \end{aligned}

Stan 代码如下

data {
int<lower=0> T;   // # time points (equally spaced)
vector[T] y;      // mean corrected return at time t
}
parameters {
real mu;                     // mean log volatility
real<lower=-1, upper=1> phi; // persistence of volatility
real<lower=0> sigma;         // white noise shock scale
vector[T] h_std;             // std log volatility time t
}
transformed parameters {
vector[T] h = h_std * sigma;  // now h ~ normal(0, sigma)
h[1] /= sqrt(1 - phi * phi);  // rescale h[1]
h += mu;
for (t in 2:T) {
h[t] += phi * (h[t - 1] - mu);
}
}
model {
phi ~ uniform(-1, 1);
sigma ~ cauchy(0, 5);
mu ~ cauchy(0, 10);

h_std ~ std_normal();
y ~ normal(0, exp(h / 2));
}

library(cmdstanr)
# 编译模型
mod_volatility_normal <- cmdstan_model(
stan_file = "code/stochastic_volatility_models.stan",
compile = TRUE, cpp_options = list(stan_threads = TRUE)
)
# 准备数据
mdata = list(T = 1274, y = as.vector(meituan_log_return))
# 拟合模型
fit_volatility_normal <- mod_volatility_normal$sample( data = mdata, chains = 2, parallel_chains = 2, iter_warmup = 1000, iter_sampling = 1000, threads_per_chain = 2, seed = 20232023, show_messages = FALSE, refresh = 0 ) # 输出结果 fit_volatility_normal$summary(c("mu", "phi", "sigma", "lp__"))
#> # A tibble: 4 × 10
#>   variable     mean   median      sd     mad       q5      q95  rhat ess_bulk
#>   <chr>       <dbl>    <dbl>   <dbl>   <dbl>    <dbl>    <dbl> <dbl>    <dbl>
#> 1 mu         -6.86    -6.86   0.117   0.115    -7.05    -6.66   1.00    1322.
#> 2 phi         0.915    0.922  0.0348  0.0300    0.849    0.960  1.00     374.
#> 3 sigma       0.294    0.286  0.0627  0.0608    0.207    0.405  1.01     366.
#> 4 lp__     3088.    3087.    29.2    28.8    3041.    3136.     1.00     544.
#> # ℹ 1 more variable: ess_tail <dbl>

11.1.2 fGarch 包

《金融时间序列分析讲义》两个波动率建模方法

• 自回归条件异方差模型（Autoregressive Conditional Heteroskedasticity，简称 ARCH）。
• 广义自回归条件异方差模型 （Generalized Autoregressive Conditional Heteroskedasticity，简称 GARCH ）

acf((meituan_log_return - mean(meituan_log_return))^2, main = "")
pacf((meituan_log_return - mean(meituan_log_return))^2, main = "")

\begin{aligned} r_t &= \mu + a_t, \quad a_t = \sigma_t \epsilon_t, \quad \epsilon_t \sim \mathcal{N}(0,1) \\ \sigma_t^2 &= \alpha_0 + \alpha_1 a_{t-1}^2 + \alpha_2 a_{t-2}^2. \end{aligned}

library(fGarch)
meituan_garch1 <- garchFit(
formula = ~ 1 + garch(2, 0),
data = meituan_log_return, trace = FALSE, cond.dist = "std"
)
summary(meituan_garch1)
#>
#> Title:
#>  GARCH Modelling
#>
#> Call:
#>  garchFit(formula = ~1 + garch(2, 0), data = meituan_log_return,
#>     cond.dist = "std", trace = FALSE)
#>
#> Mean and Variance Equation:
#>  data ~ 1 + garch(2, 0)
#> <environment: 0x559a16a757c8>
#>  [data = meituan_log_return]
#>
#> Conditional Distribution:
#>  std
#>
#> Coefficient(s):
#>         mu       omega      alpha1      alpha2       shape
#> 0.00025765  0.00107291  0.11199514  0.13829052  4.93560791
#>
#> Std. Errors:
#>  based on Hessian
#>
#> Error Analysis:
#>         Estimate  Std. Error  t value Pr(>|t|)
#> mu     2.577e-04   8.970e-04    0.287  0.77394
#> omega  1.073e-03   9.432e-05   11.375  < 2e-16 ***
#> alpha1 1.120e-01   4.292e-02    2.609  0.00907 **
#> alpha2 1.383e-01   4.725e-02    2.927  0.00343 **
#> shape  4.936e+00   7.008e-01    7.043 1.88e-12 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Log Likelihood:
#>  2459.345    normalized:  1.930412
#>
#> Description:
#>  Wed Jun 19 23:52:12 2024 by user:
#>
#>
#> Standardised Residuals Tests:
#>                                  Statistic      p-Value
#>  Jarque-Bera Test   R    Chi^2  260.648221 0.000000e+00
#>  Shapiro-Wilk Test  R    W        0.975911 9.515068e-14
#>  Ljung-Box Test     R    Q(10)   21.212805 1.965757e-02
#>  Ljung-Box Test     R    Q(15)   24.773613 5.306761e-02
#>  Ljung-Box Test     R    Q(20)   33.252172 3.165156e-02
#>  Ljung-Box Test     R^2  Q(10)   17.362615 6.671552e-02
#>  Ljung-Box Test     R^2  Q(15)   29.329025 1.458471e-02
#>  Ljung-Box Test     R^2  Q(20)   53.703424 6.400354e-05
#>  LM Arch Test       R    TR^2    19.254099 8.257804e-02
#>
#> Information Criterion Statistics:
#>       AIC       BIC       SIC      HQIC
#> -3.852975 -3.832764 -3.853006 -3.845384

\begin{aligned} r_t &= -5.665 \times 10^{-5} + a_t, \quad a_t = \sigma_t \epsilon_t, \quad \epsilon_t \sim \mathcal{N}(0,1) \\ \sigma_t^2 &= 1.070 \times 10^{-3} + 0.1156 a_{t-1}^2 + 0.1438a_{t-2}^2. \end{aligned}

\begin{aligned} r_t &= \mu + a_t, \quad a_t = \sigma_t \epsilon_t, \quad \epsilon_t \sim \mathcal{N}(0,1) \\ \sigma_t^2 &= \alpha_0 + \alpha_1 a_{t-1}^2 + \beta_1 \sigma_{t-1}^2. \end{aligned}

meituan_garch2 <- garchFit(
formula = ~ 1 + garch(1, 1),
data = meituan_log_return, trace = FALSE, cond.dist = "std"
)
summary(meituan_garch2)
#>
#> Title:
#>  GARCH Modelling
#>
#> Call:
#>  garchFit(formula = ~1 + garch(1, 1), data = meituan_log_return,
#>     cond.dist = "std", trace = FALSE)
#>
#> Mean and Variance Equation:
#>  data ~ 1 + garch(1, 1)
#> <environment: 0x559a0b3a8c60>
#>  [data = meituan_log_return]
#>
#> Conditional Distribution:
#>  std
#>
#> Coefficient(s):
#>         mu       omega      alpha1       beta1       shape
#> 2.8294e-04  3.4453e-05  5.9798e-02  9.1678e-01  5.4352e+00
#>
#> Std. Errors:
#>  based on Hessian
#>
#> Error Analysis:
#>         Estimate  Std. Error  t value Pr(>|t|)
#> mu     2.829e-04   8.702e-04    0.325  0.74507
#> omega  3.445e-05   1.937e-05    1.779  0.07525 .
#> alpha1 5.980e-02   1.855e-02    3.224  0.00127 **
#> beta1  9.168e-01   2.784e-02   32.933  < 2e-16 ***
#> shape  5.435e+00   8.137e-01    6.680 2.39e-11 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Log Likelihood:
#>  2478.473    normalized:  1.945427
#>
#> Description:
#>  Wed Jun 19 23:52:12 2024 by user:
#>
#>
#> Standardised Residuals Tests:
#>                                   Statistic      p-Value
#>  Jarque-Bera Test   R    Chi^2  226.7201557 0.000000e+00
#>  Shapiro-Wilk Test  R    W        0.9781165 5.582407e-13
#>  Ljung-Box Test     R    Q(10)   16.0489321 9.824020e-02
#>  Ljung-Box Test     R    Q(15)   19.6491167 1.858101e-01
#>  Ljung-Box Test     R    Q(20)   27.2460656 1.284820e-01
#>  Ljung-Box Test     R^2  Q(10)    7.8054889 6.478299e-01
#>  Ljung-Box Test     R^2  Q(15)    9.8336965 8.300684e-01
#>  Ljung-Box Test     R^2  Q(20)   24.7640835 2.106061e-01
#>  LM Arch Test       R    TR^2     9.6000051 6.510060e-01
#>
#> Information Criterion Statistics:
#>       AIC       BIC       SIC      HQIC
#> -3.883004 -3.862792 -3.883034 -3.875413

11.2 贝叶斯可加模型

$y = at + b + c \sin(\frac{t}{12} \times 2\pi) + d \cos(\frac{t}{12} \times 2\pi) + \epsilon$

air_passengers_df <- data.frame(y = as.vector(AirPassengers), t = 1:144)
fit_lm1 <- lm(y ~ t + sin(t / 12 * 2 * pi) + cos(t / 12 * 2 * pi), data = air_passengers_df)
fit_lm2 <- update(fit_lm1, . ~ . +
sin(t / 12 * 2 * 2 * pi) + cos(t / 12 * 2 * 2 * pi), data = air_passengers_df
)
fit_lm3 <- update(fit_lm2, . ~ . +
sin(t / 12 * 3 * 2 * pi) + cos(t / 12 * 3 * 2 * pi), data = air_passengers_df
)
plot(y ~ t, air_passengers_df, type = "l")
lines(x = air_passengers_df$t, y = fit_lm1$fitted.values, col = "red")
lines(x = air_passengers_df$t, y = fit_lm2$fitted.values, col = "green")
lines(x = air_passengers_df$t, y = fit_lm3$fitted.values, col = "orange")

11.2.1 Stan 框架

prophet 包是如何同时处理这些情况，是否可以在 cmdstanr 包中实现，是否可以在 mgcv 和 INLA 中实现？

library(cmdstanr)

11.2.2 INLA 框架

library(ggfortify)
autoplot(log(AirPassengers)) +
theme_classic() +
labs(x = "年月", y = "对数值")
autoplot(diff(log(AirPassengers))) +
theme_classic() +
labs(x = "年月", y = "差分对数值")

library(INLA)
inla.setOption(short.summary = TRUE)
air_passengers_df <- data.frame(
y = as.vector(AirPassengers),
year = as.factor(rep(1949:1960, each = 12)),
month = as.factor(rep(1:12, times = 12)),
ID = 1:length(AirPassengers)
)
mod_inla_rw1 <- inla(
formula = log(y) ~ year + f(ID, model = "rw1"),
family = "gaussian", data = air_passengers_df,
control.family = list(hyper = list(prec = list(param = c(1, 0.2)))),
control.predictor = list(compute = TRUE)
)
summary(mod_inla_rw1)
#> Fixed effects:
#>              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
#> (Intercept) 5.159 0.252      4.666    5.158      5.657 5.158   0
#> year1950    0.050 0.134     -0.215    0.050      0.313 0.050   0
#> year1951    0.174 0.190     -0.201    0.174      0.546 0.174   0
#> year1952    0.262 0.233     -0.198    0.262      0.717 0.262   0
#> year1953    0.330 0.269     -0.201    0.331      0.857 0.331   0
#> year1954    0.357 0.301     -0.236    0.358      0.945 0.358   0
#> year1955    0.442 0.329     -0.208    0.443      1.086 0.443   0
#> year1956    0.510 0.356     -0.193    0.511      1.206 0.511   0
#> year1957    0.567 0.381     -0.184    0.568      1.311 0.568   0
#> year1958    0.576 0.403     -0.220    0.577      1.365 0.577   0
#> year1959    0.647 0.425     -0.191    0.648      1.478 0.648   0
#> year1960    0.683 0.445     -0.195    0.684      1.555 0.684   0
#>
#> Model hyperparameters:
#>                                           mean    sd 0.025quant 0.5quant
#> Precision for the Gaussian observations 101.38 18.14      69.78    99.99
#> Precision for ID                        155.14 40.69      92.69   149.32
#>                                         0.975quant   mode
#> Precision for the Gaussian observations     140.94  97.59
#> Precision for ID                            251.66 137.50
#>
#>  is computed

mod_inla_sea <- inla(
formula = log(y) ~ year + f(ID, model = "seasonal", season.length = 12),
family = "gaussian", data = air_passengers_df,
control.family = list(hyper = list(prec = list(param = c(1, 0.2)))),
control.predictor = list(compute = TRUE)
)
summary(mod_inla_sea)
#> Fixed effects:
#>              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
#> (Intercept) 4.836 0.020      4.797    4.836      4.875 4.836   0
#> year1950    0.095 0.028      0.039    0.095      0.150 0.095   0
#> year1951    0.295 0.028      0.240    0.295      0.351 0.295   0
#> year1952    0.441 0.028      0.386    0.441      0.496 0.441   0
#> year1953    0.573 0.028      0.517    0.573      0.628 0.573   0
#> year1954    0.630 0.028      0.575    0.630      0.686 0.630   0
#> year1955    0.803 0.028      0.748    0.803      0.859 0.803   0
#> year1956    0.948 0.028      0.893    0.948      1.004 0.948   0
#> year1957    1.062 0.028      1.007    1.062      1.118 1.062   0
#> year1958    1.094 0.028      1.039    1.094      1.150 1.094   0
#> year1959    1.212 0.028      1.157    1.212      1.268 1.212   0
#> year1960    1.318 0.028      1.263    1.318      1.373 1.318   0
#>
#> Model hyperparameters:
#>                                             mean       sd 0.025quant 0.5quant
#> Precision for the Gaussian observations   213.03    27.48     163.46   211.48
#> Precision for ID                        42079.80 27524.61   10127.06 35355.57
#>                                         0.975quant     mode
#> Precision for the Gaussian observations     271.44   208.83
#> Precision for ID                         113478.54 24333.72
#>
#>  is computed

mod_inla_rw1_fitted <- data.frame(
ID = 1:length(AirPassengers),
y = as.vector(log(AirPassengers)),
mean = mod_inla_rw1$summary.fitted.values$mean,
0.025quant = mod_inla_rw1$summary.fitted.values$0.025quant,
0.975quant = mod_inla_rw1$summary.fitted.values$0.975quant,
check.names = FALSE
)
mod_inla_sea_fitted <- data.frame(
ID = 1:length(AirPassengers),
y = as.vector(log(AirPassengers)),
mean = mod_inla_sea$summary.fitted.values$mean,
0.025quant = mod_inla_sea$summary.fitted.values$0.025quant,
0.975quant = mod_inla_sea$summary.fitted.values$0.975quant,
check.names = FALSE
)
ggplot(data = mod_inla_rw1_fitted, aes(ID)) +
geom_ribbon(aes(ymin = 0.025quant, ymax = 0.975quant), fill = "gray") +
geom_line(aes(y = y)) +
geom_line(aes(y = mean), color = "red") +
theme_classic() +
labs(x = "序号", y = "对数值")
ggplot(data = mod_inla_sea_fitted, aes(ID)) +
geom_ribbon(aes(ymin = 0.025quant, ymax = 0.975quant), fill = "gray") +
geom_line(aes(y = y)) +
geom_line(aes(y = mean), color = "red") +
theme_classic() +
labs(x = "序号", y = "对数值")

11.3 一些非参数模型

11.3.1 mgcv 包

mgcv 是 R 软件内置的推荐组件，由 Simon Wood 开发和维护，历经多年，成熟稳定。函数 bam() 相比于函数 gam() 的优势是可以处理大规模的时间序列数据。对于时间序列数据预测，数万和百万级观测值都可以

air_passengers_tbl <- data.frame(
y = as.vector(AirPassengers),
year = rep(1949:1960, each = 12),
month = rep(1:12, times = 12)
)
mod1 <- gam(y ~ s(year) + s(month, bs = "cr", k = 12),
data = air_passengers_tbl, family = gaussian
)
summary(mod1)
#>
#> Family: gaussian
#>
#> Formula:
#> y ~ s(year) + s(month, bs = "cr", k = 12)
#>
#> Parametric coefficients:
#>             Estimate Std. Error t value Pr(>|t|)
#> (Intercept)  280.299      1.957   143.2   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Approximate significance of smooth terms:
#>            edf Ref.df      F p-value
#> s(year)  6.102  7.265 441.40  <2e-16 ***
#> s(month) 8.796 10.097  38.25  <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> R-sq.(adj) =  0.962   Deviance explained = 96.6%
#> GCV =  619.9  Scale est. = 551.47    n = 144

layout(matrix(1:2, nrow = 1))
plot(mod1, shade = TRUE)

air_passengers_ts <- ts(mod1$fitted.values, start = c(1949, 1), frequency = 12) plot(AirPassengers) lines(air_passengers_ts, col = "red") 整体上，乘客数逐年呈线性增长，每年不同月份呈现波动，淡季和旺季出行的流量有很大差异，近年来，这种差异的波动在扩大。为了刻画这种情况，考虑年度趋势和月度波动的交互作用。 mod2 <- gam(y ~ s(year, month), data = air_passengers_tbl, family = gaussian) summary(mod2) #> #> Family: gaussian #> Link function: identity #> #> Formula: #> y ~ s(year, month) #> #> Parametric coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 280.299 1.059 264.7 <2e-16 *** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> Approximate significance of smooth terms: #> edf Ref.df F p-value #> s(year,month) 28.21 28.96 435.3 <2e-16 *** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> R-sq.(adj) = 0.989 Deviance explained = 99.1% #> GCV = 202.62 Scale est. = 161.52 n = 144 可以看到，调整的 $$R^2$$ 明显增加，拟合效果更好，各年各月份的乘客数变化，见下图。 op <- par(mar = c(4, 4, 2, 0)) plot(mod2) on.exit(par(op), add = TRUE)  上图是轮廓图，下面用透视图展示趋势拟合的效果。 op <- par(mar = c(0, 1.5, 0, 0)) vis.gam(mod2, theta = -35, phi = 20, ticktype = "detailed", expand = .65, zlab = "") on.exit(par(op), add = TRUE)  最后，在原始数据的基础上，添加拟合数据，得到如下拟合趋势图，与前面的拟合图比较，可以看出效果提升很明显。 air_passengers_ts <- ts(mod2$fitted.values, start = c(1949, 1), frequency = 12)
plot(AirPassengers)
lines(air_passengers_ts, col = "red")

11.3.2 nnet 包

# 准备数据
air_passengers <- as.matrix(embed(AirPassengers, 4))
colnames(air_passengers) <- c("y", "x3", "x2", "x1")
data_size <- nrow(air_passengers)
# 拆分数据集
train_size <- floor(data_size * 0.67)
train_data <- air_passengers[1:train_size, ]
test_data <- air_passengers[train_size:data_size, ]

# 随机数种子对结果的影响非常大 试试 set.seed(20232023)
set.seed(20222022)
# 单隐藏层 8 个神经元
mod_nnet <- nnet::nnet(
y ~ x1 + x2 + x3,
data = air_passengers, # 数据集
subset = 1:train_size, # 训练数据的指标向量
linout = TRUE, size = 4, rang = 0.1,
decay = 5e-4, maxit = 400, trace = FALSE
)
# 预测
train_pred <- predict(mod_nnet, newdata = air_passengers[1:train_size,], type = "raw")
# 训练集 RMSE
sqrt(mean((air_passengers[1:train_size, "y"] - train_pred )^2))
#> [1] 21.59392
# 预测
test_pred <- predict(mod_nnet, newdata = air_passengers[-(1:train_size),], type = "raw")
# 测试集 RMSE
sqrt(mean((air_passengers[-(1:train_size), "y"] - test_pred)^2))
#> [1] 53.79107

train_pred_ts <- ts(data = train_pred, start = c(1949, 3), frequency = 12)
test_pred_ts <- ts(data = test_pred, start = c(1957, 1), frequency = 12)
plot(AirPassengers)
lines(train_pred_ts, col = "red")
lines(test_pred_ts, col = "green")

11.3.3 keras3 包

keras3 包通过 reticulate 包引入 Keras 3 框架，这个框架支持 TensorFlowPyTorch 等多个后端，目前，keras3 包通过 tensorflow 包仅支持 TensorFlow 后端。

library(keras3)
set_random_seed(20222022)
# 模型结构
mod_mlp <- keras_model_sequential(shape = c(3)) |>
layer_dense(units = 12, activation = "relu") |>
layer_dense(units = 8, activation = "relu") |>
layer_dense(units = 1)
# 训练目标
compile(mod_mlp,
loss = "mse", # 损失函数
metrics = "mae" # 监控度量
)
# 模型概览
summary(mod_mlp)
#> Model: "sequential"
#> ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┓
#> ┃ Layer (type)                      ┃ Output Shape             ┃       Param # ┃
#> ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━┩
#> │ dense (Dense)                     │ (None, 12)               │            48 │
#> ├───────────────────────────────────┼──────────────────────────┼───────────────┤
#> │ dense_1 (Dense)                   │ (None, 8)                │           104 │
#> ├───────────────────────────────────┼──────────────────────────┼───────────────┤
#> │ dense_2 (Dense)                   │ (None, 1)                │             9 │
#> └───────────────────────────────────┴──────────────────────────┴───────────────┘
#>  Total params: 161 (644.00 B)
#>  Trainable params: 161 (644.00 B)
#>  Non-trainable params: 0 (0.00 B)

# 拟合模型
fit(mod_mlp,
x = train_data[, c("x1", "x2", "x3")],
y = train_data[, "y"],
epochs = 200,
batch_size = 10, # 每次更新梯度所用的样本量
validation_split = 0.2, # 从训练数据中拆分一部分用作验证集
verbose = 0 # 不显示训练进度
)
# 将测试数据代入模型，计算损失函数和监控度量
evaluate(mod_mlp, test_data[, c("x1", "x2", "x3")], test_data[, "y"])
#> 2/2 - 0s - 8ms/step - loss: 2529.3757 - mae: 41.1564
#> $loss #> [1] 2529.376 #> #>$mae
#> [1] 41.15638
# 测试集上的预测
mlp_test_pred <- predict(mod_mlp, test_data[, c("x1", "x2", "x3")]) 
#> 2/2 - 0s - 28ms/step
mlp_train_pred <- predict(mod_mlp, train_data[, c("x1", "x2", "x3")]) 
#> 3/3 - 0s - 4ms/step
sqrt(mean((test_data[, "y"] - mlp_test_pred)^2)) # 计算均方根误差
#> [1] 50.2929

mlp_train_pred_ts <- ts(data = mlp_train_pred, start = c(1949, 3), frequency = 12)
mlp_test_pred_ts <- ts(data = mlp_test_pred, start = c(1957, 1), frequency = 12)
plot(AirPassengers)
lines(mlp_train_pred_ts, col = "red")
lines(mlp_test_pred_ts, col = "green")

11.4 习题

1. 基于 R 软件内置的数据集 sunspotssunspot.month 比较 INLA 和 mgcv 框架的预测效果。

代码
sunspots_tbl <- broom::tidy(sunspots)
sunspots_month_tbl <- broom::tidy(sunspot.month)
ggplot() +
geom_line(data = sunspots_month_tbl, aes(x = index, y = value), color = "red") +
geom_line(data = sunspots_tbl, aes(x = index, y = value)) +
theme_bw() +
labs(x = "年月", y = "数量")

图中黑线和红线分别表示 1749-1983 年、1984-2014 年每月太阳黑子数量。