\[\rm{NIKKEI}\sim\rm{Normal}(\beta_0 + \beta_1 \cdot \rm{USDJPY},\sigma)\]

# 数値はいずれも前月比(%)
head(df)
        Date NIKKEI USDJPY
2 2017-02-01  -0.03  -1.45
3 2017-03-01   0.79  -0.05
4 2017-04-01  -3.12  -2.59
5 2017-05-01   5.29   1.94
6 2017-06-01   1.62  -1.19
7 2017-07-01   0.00   1.36
tail(df)
         Date NIKKEI USDJPY
32 2019-08-01  -4.46  -1.81
33 2019-09-01   4.63   1.06
34 2019-10-01   2.84   0.65
35 2019-11-01   4.87   0.69
36 2019-12-01   1.64   0.32
37 2020-01-01  -0.61  -0.55
apply(df[, -1], 2, adf.test)
$NIKKEI

    Augmented Dickey-Fuller Test

data:  newX[, i]
Dickey-Fuller = -3.545, Lag order = 3, p-value = 0.05218
alternative hypothesis: stationary


$USDJPY

    Augmented Dickey-Fuller Test

data:  newX[, i]
Dickey-Fuller = -3.7536, Lag order = 3, p-value = 0.03533
alternative hypothesis: stationary
# 最尤推定
summary(lm(NIKKEI ~ USDJPY), confint = T, ci.width = 0.95)

Call:
lm(formula = NIKKEI ~ USDJPY)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.7426 -1.0441  0.2489  1.1406  6.1535 

Coefficients:
            Estimate Std. Error t value  Pr(>|t|)    
(Intercept)   0.8275     0.4212   1.965    0.0577 .  
USDJPY        1.4797     0.2977   4.970 0.0000188 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.514 on 34 degrees of freedom
Multiple R-squared:  0.4208,    Adjusted R-squared:  0.4038 
F-statistic:  24.7 on 1 and 34 DF,  p-value: 0.00001878
Gaussian <- "
  data{
    int N;
    vector[N] NIKKEI;
    vector[N] USDJPY;
  }
  parameters{
    real beta0;
    real beta1;
    real <lower = 0> sigma;
  }
  model{
    for(i in 1:N)
      NIKKEI[i] ~ normal(beta0 + beta1*USDJPY[i], sigma);
  }
  generated quantities{
    vector[N] pred_NIKKEI;
    real log_lik[N];
    for (i in 1:N){
      pred_NIKKEI[i] = normal_rng(beta0 + beta1*USDJPY[i], sigma);
      log_lik[i] = normal_lpdf(NIKKEI[i] | beta0 + beta1*USDJPY[i], sigma);
    }
  }
"
datalist <- list(N = N, NIKKEI = NIKKEI, USDJPY = USDJPY)
iter <- 1400
warmup <- 400
chains <- 3
fit <- stan(model_code = Gaussian, data = datalist, iter = iter, warmup = warmup, thin = 1, chains = chains)
summary(fit)$summary[c("beta0", "beta1", "sigma"), ]
           mean     se_mean        sd        2.5%       25%      50%      75%    97.5%    n_eff      Rhat
beta0 0.8157041 0.008776048 0.4329639 -0.05810871 0.5339279 0.820452 1.106948 1.652507 2433.915 1.0010500
beta1 1.4790183 0.005919812 0.3101240  0.86576085 1.2765948 1.474421 1.680555 2.110719 2744.447 0.9997409
sigma 2.6107961 0.006883586 0.3346165  2.06597912 2.3633772 2.578157 2.818477 3.370859 2363.008 1.0005781
traceplot(fit) + theme(axis.text.x = element_text(size = 5), axis.text.y = element_text(size = 5), strip.text.x = element_text(size = 5), legend.title = element_text(size = 5), legend.text = element_text(size = 5))

# EAP:事後期待値
df_result <- rstan::extract(fit)$pred_NIKKEI %>% data.frame() %>% gather() %>% dplyr::mutate(id = rep(c(1:N), each = (iter - warmup) * chains)) %>% group_by(id) %>% dplyr::summarize(pred_EAP = mean(value), pred_lower = quantile(value, 0.025), pred_upper = quantile(value, 0.975)) %>% dplyr::ungroup() %>% cbind(data.frame(NIKKEI, USDJPY))
head(df_result)
  id   pred_EAP pred_lower pred_upper NIKKEI USDJPY
1  1 -1.3095587  -6.622188   3.992890  -0.03  -1.45
2  2  0.6987707  -4.569086   5.862750   0.79  -0.05
3  3 -3.1180798  -8.678892   2.415784  -3.12  -2.59
4  4  3.6821586  -1.816969   9.053668   5.29   1.94
5  5 -0.8368877  -6.304224   4.354736   1.62  -1.19
6  6  2.8763132  -2.601950   8.016950   0.00   1.36
tail(df_result)
   id    pred_EAP pred_lower pred_upper NIKKEI USDJPY
31 31 -1.80668264  -7.387550   3.851648  -4.46  -1.81
32 32  2.42867352  -2.825009   7.576436   4.63   1.06
33 33  1.80444039  -3.556048   7.058846   2.84   0.65
34 34  1.92284211  -3.245637   7.230945   4.87   0.69
35 35  1.25941029  -4.040424   6.449707   1.64   0.32
36 36  0.04027954  -5.304679   5.447900  -0.61  -0.55