forked from jrnold/bugs-examples-in-stan
-
Notifications
You must be signed in to change notification settings - Fork 1
/
reagan.Rmd
92 lines (78 loc) · 2.79 KB
/
reagan.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
# Reagan: linear regression with AR(1) disturbances {#reagan}
```{r reagan_setup,message=FALSE,cache=FALSE}
library("tidyverse")
library("rstan")
```
Ninety-six monthly observations on presidential job approval ratings for Ronald Reagan are modeled via linear regression, with a correction for first-order serial correlation among the disturbances.[^reagan]
Note the marginal model for the first observation, and the conditioning on the lagged observation for months 2 through 96.
A uniform prior over the stationary (-1,1) interval is employed for the residual AR(1) parameter.
$$
\begin{aligned}[t]
y_i &= \mu_i + \epsilon_i + \theta \epsilon_{i - 1} ,\\
\mu_i &= \alpha + x_i' \beta , \\
\epsilon_i &\sim \mathsf{Normal}(0, \sigma^2) ,
\end{aligned}
$$
for $i \in 1, \dots, N$.
Weakly informative priors for each parameter are used,
$$
\begin{aligned}[t]
\alpha &\sim \mathsf{Normal}(0, 10), \\
\beta_k &\sim \mathsf{Normal}(0, 2.5), & k \in 1, \dots, K, \\
\sigma &\sim \mathsf{HalfCauchy}(0, 5), \\
\theta &= 2 \theta^{*} - 1 , \\
\theta^{*} &\sim \mathsf{Beta}(1, 1) .
\end{aligned}
$$
```{r ReaganApproval}
data("ReaganApproval", package = "bayesjackman")
ReaganApproval
```
```{r reagan_data}
reagan_data <- within(list(), {
y <- ReaganApproval$app
N <- length(y)
X <- model.matrix(~ 0 + infl + unemp, data = ReaganApproval) %>% scale()
K <- ncol(X)
alpha_loc <- 0
alpha_scale <- 10
beta_loc <- rep(0, K)
beta_scale <- rep(2.5 * sd(y), K)
sigma_scale <- 5 * sd(y)
theta_a <- 1
theta_b <- 1
})
```
```{r mod_regar1,cache.extra=tools::md5sum("stan/regar1.stan")}
mod_regar1 <- stan_model("stan/regar1.stan")
```
```{r}
mod_regar1
```
```{r reagan_fit,results='hide'}
reagan_fit <- sampling(mod_regar1, data = reagan_data)
```
```{r reagan_fit_summary}
summary(reagan_fit, par = c("alpha", "beta", "theta", "sigma"))$summary
```
## Cochrane-Orcutt/Prais-Winsten
An AR(1) error model can also be estimated [Prais-Winsten](https://en.wikipedia.org/wiki/Prais%E2%80%93Winsten_estimation) estimation:
$$
\begin{aligned}[t]
y_1 &\sim \mathsf{Normal}\left(\alpha + x_1' \beta, \frac{\sigma ^ 2}{1 - \theta ^ 2} \right), \\
y_i &\sim \mathsf{Normal}\left(\theta y_{i - 1} + \alpha (1 - \theta) + \beta (X_i - \theta X_{i - 1}), \sigma ^ 2 \right) & i = 2, \dots, N
\end{aligned}
$$
```{r mod_pw,cache.extra=tools::md5sum("stan/pw.stan")}
mod_pw <- stan_model("stan/pw.stan")
```
```{r}
mod_pw
```
```{r reagan_fit2,results='hide'}
reagan_fit2 <- sampling(mod_pw, data = reagan_data)
```
```{r reagan_fit_summary2}
summary(reagan_fit2, par = c("alpha", "beta", "theta", "sigma"))$summary
```
[^reagan]: Example derived from Simon Jackman, "Reagan: linear regression with AR(1) disturbances," *BUGS Examples,* 2007-07-24, [URL](https://web-beta.archive.org/web/20070724034151/http://jackman.stanford.edu:80/mcmc/reagan.odc).