Bayesian Regression in Rstan

Loading Data

opts_chunk$set(warning = FALSE, error = FALSE, message = FALSE)
require(plm)
require(gdata)
require(DMwR)
require(rstan)
setwd("/Users/Acton/Desktop/adaFdata")
CPI = read.xls(skip = 50, xls = "CPIAUCSL (3).xls")
DGS10 = read.xls(skip = 14, xls = "DGS10.xls")  #10-Y Treasury Constant Maturity Rate
GDP = read.xls(skip = 17, xls = "GDPC1.xls")
UNRATE = read.xls(skip = 18, xls = "UNRATE.xls")
SP500 = read.xls(skip = 17, xls = "SP500 (1).xls")
data = cbind(SP500, CPI, DGS10, GDP, UNRATE)[, -c(3, 5, 7, 9)]
colnames(data) = c("Date", "SP500", "CPI", "DGS10", "GDP", "UNRATE")

Data Tranformation

Transfer data to get rid and serial correlation in time series data and scale it

logdata = log(data[, -1])
scaledata = as.data.frame(scale(logdata))
dfdata = as.data.frame(diff(as.matrix(scaledata), lag = 1, differenc = 1))

Fit Fix Effect Model

To fit a fix effect model, I will take advantage of existing r function. Here I treat the data as time series cross sectional data in some fashion. And only account for individual effect(which is the year effect in my data,and the time effect is season effect)

pdata = pdata.frame(dfdata[1:92, ], index = 23, drop.index = T, row.names = T)
pm = plm(SP500 ~ CPI + DGS10 + GDP + UNRATE - 1, data = pdata, effect = c("individual"), 
    model = "within")
fixef(pm, effect = "individual")
##         1         2         3         4         5         6         7 
## -0.028039  0.014589 -0.027216 -0.039652 -0.079688  0.100831  0.007124 
##         8         9        10        11        12        13        14 
##  0.079163  0.037236 -0.052262 -0.104610 -0.103668 -0.175167  0.062770 
##        15        16        17        18        19        20        21 
## -0.060521 -0.048316  0.008024 -0.071781 -0.205388  0.089875  0.034232 
##        22        23 
## -0.005846  0.028824

Fit Bayesian Model Accounting for Year Effect

I = matrix(1, 4, 1)
I1 = diag(1, 23, 23)
Iyear1 = cbind(I1 %x% I, matrix(0, 92, 1))
Iyear2 = cbind(matrix(0, 3, 23), matrix(1, 3, 1))
Iyear = rbind(Iyear1, Iyear2)
dfdataexp <- cbind(dfdata, Iyear)
BM <- "\ndata {\nint<lower=0> N;\nint<lower=0> J;\nint<lower=0> J1;\nreal y[N];\nmatrix[N,J] x;\nmatrix[N,J1] year;\n}\nparameters {\nvector[J] w;\nvector[J1] wyear;\nreal<lower=0> S;\n}\ntransformed parameters{\nvector[N] U; \nU<-x * w + year * wyear;\n}\n\nmodel {\nfor(j1 in 1:J1)\nyear[j1] ~ normal(0,1);\nfor(j in 1:J)\nw[j] ~ normal(0,10);\nS ~ inv_gamma(1,2);\nfor(n in 1:N)\ny[n] ~ normal(U[n],S);\n}\n"
Bdata <- list(J = 4, N = 95, y = dfdataexp[, 1], x = dfdataexp[, c(2:5)], year = dfdataexp[, 
    -c(1:5)], J1 = 24)
SP_fit <- stan(model_code = BM, data = Bdata, iter = 600, chains = 4)
## 
## TRANSLATING MODEL 'BM' FROM Stan CODE TO C++ CODE NOW.
## COMPILING THE C++ CODE FOR MODEL 'BM' NOW.
## SAMPLING FOR MODEL 'BM' NOW (CHAIN 1).
## 
Iteration:   1 / 600 [  0%]  (Warmup)
Iteration:  60 / 600 [ 10%]  (Warmup)
Iteration: 120 / 600 [ 20%]  (Warmup)
Iteration: 180 / 600 [ 30%]  (Warmup)
Iteration: 240 / 600 [ 40%]  (Warmup)
Iteration: 300 / 600 [ 50%]  (Warmup)
Iteration: 360 / 600 [ 60%]  (Sampling)
Iteration: 420 / 600 [ 70%]  (Sampling)
Iteration: 480 / 600 [ 80%]  (Sampling)
Iteration: 540 / 600 [ 90%]  (Sampling)
Iteration: 600 / 600 [100%]  (Sampling)
## Elapsed Time: 3.79067 seconds (Warm-up)
##               3.07445 seconds (Sampling)
##               6.86512 seconds (Total)
## 
## SAMPLING FOR MODEL 'BM' NOW (CHAIN 2).
## 
Iteration:   1 / 600 [  0%]  (Warmup)
Iteration:  60 / 600 [ 10%]  (Warmup)
Iteration: 120 / 600 [ 20%]  (Warmup)
Iteration: 180 / 600 [ 30%]  (Warmup)
Iteration: 240 / 600 [ 40%]  (Warmup)
Iteration: 300 / 600 [ 50%]  (Warmup)
Iteration: 360 / 600 [ 60%]  (Sampling)
Iteration: 420 / 600 [ 70%]  (Sampling)
Iteration: 480 / 600 [ 80%]  (Sampling)
Iteration: 540 / 600 [ 90%]  (Sampling)
Iteration: 600 / 600 [100%]  (Sampling)
## Elapsed Time: 3.97697 seconds (Warm-up)
##               3.93525 seconds (Sampling)
##               7.91222 seconds (Total)
## 
## SAMPLING FOR MODEL 'BM' NOW (CHAIN 3).
## 
Iteration:   1 / 600 [  0%]  (Warmup)
Iteration:  60 / 600 [ 10%]  (Warmup)
Iteration: 120 / 600 [ 20%]  (Warmup)
Iteration: 180 / 600 [ 30%]  (Warmup)
Iteration: 240 / 600 [ 40%]  (Warmup)
Iteration: 300 / 600 [ 50%]  (Warmup)
Iteration: 360 / 600 [ 60%]  (Sampling)
Iteration: 420 / 600 [ 70%]  (Sampling)
Iteration: 480 / 600 [ 80%]  (Sampling)
Iteration: 540 / 600 [ 90%]  (Sampling)
Iteration: 600 / 600 [100%]  (Sampling)
## Elapsed Time: 3.81817 seconds (Warm-up)
##               4.08911 seconds (Sampling)
##               7.90728 seconds (Total)
## 
## SAMPLING FOR MODEL 'BM' NOW (CHAIN 4).
## 
Iteration:   1 / 600 [  0%]  (Warmup)
Iteration:  60 / 600 [ 10%]  (Warmup)
Iteration: 120 / 600 [ 20%]  (Warmup)
Iteration: 180 / 600 [ 30%]  (Warmup)
Iteration: 240 / 600 [ 40%]  (Warmup)
Iteration: 300 / 600 [ 50%]  (Warmup)
Iteration: 360 / 600 [ 60%]  (Sampling)
Iteration: 420 / 600 [ 70%]  (Sampling)
Iteration: 480 / 600 [ 80%]  (Sampling)
Iteration: 540 / 600 [ 90%]  (Sampling)
Iteration: 600 / 600 [100%]  (Sampling)
## Elapsed Time: 4.22515 seconds (Warm-up)
##               4.53629 seconds (Sampling)
##               8.76143 seconds (Total)
plot(SP_fit)

plot of chunk unnamed-chunk-9

print(SP_fit)
## Inference for Stan model: BM.
## 4 chains, each with iter=600; warmup=300; thin=1; 
## post-warmup draws per chain=300, total post-warmup draws=1200.
## 
##            mean se_mean  sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
## w[1]        0.9     0.0 0.5  -0.2   0.5   0.9   1.2   1.9   529    1
## w[2]        0.1     0.0 0.1   0.0   0.1   0.1   0.2   0.2   683    1
## w[3]        0.9     0.0 0.5  -0.2   0.6   1.0   1.3   2.0   486    1
## w[4]        0.0     0.0 0.1  -0.2  -0.1   0.0   0.1   0.3   667    1
## wyear[1]    0.0     0.0 0.1  -0.2  -0.1   0.0   0.0   0.1  1200    1
## wyear[2]    0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.1   759    1
## wyear[3]    0.0     0.0 0.1  -0.1  -0.1   0.0   0.0   0.1   740    1
## wyear[4]    0.0     0.0 0.1  -0.2  -0.1   0.0   0.0   0.1   739    1
## wyear[5]   -0.1     0.0 0.1  -0.2  -0.1  -0.1   0.0   0.1   785    1
## wyear[6]    0.1     0.0 0.1   0.0   0.1   0.1   0.1   0.2   898    1
## wyear[7]    0.0     0.0 0.1  -0.1   0.0   0.0   0.0   0.1   776    1
## wyear[8]    0.1     0.0 0.1  -0.1   0.0   0.1   0.1   0.2   567    1
## wyear[9]    0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.2   639    1
## wyear[10]   0.0     0.0 0.1  -0.2  -0.1   0.0   0.0   0.1   855    1
## wyear[11]  -0.1     0.0 0.1  -0.2  -0.2  -0.1  -0.1   0.0  1200    1
## wyear[12]  -0.1     0.0 0.1  -0.2  -0.1  -0.1  -0.1   0.0   935    1
## wyear[13]  -0.2     0.0 0.1  -0.3  -0.2  -0.2  -0.1   0.0  1200    1
## wyear[14]   0.1     0.0 0.1  -0.1   0.0   0.1   0.1   0.2  1200    1
## wyear[15]  -0.1     0.0 0.1  -0.2  -0.1  -0.1   0.0   0.1  1200    1
## wyear[16]   0.0     0.0 0.1  -0.2  -0.1   0.0   0.0   0.1   663    1
## wyear[17]   0.0     0.0 0.1  -0.1   0.0   0.0   0.0   0.1  1200    1
## wyear[18]  -0.1     0.0 0.1  -0.2  -0.1  -0.1   0.0   0.0   804    1
## wyear[19]  -0.2     0.0 0.1  -0.4  -0.3  -0.2  -0.2  -0.1   798    1
## wyear[20]   0.1     0.0 0.1   0.0   0.1   0.1   0.1   0.2  1200    1
## wyear[21]   0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.2  1200    1
## wyear[22]   0.0     0.0 0.1  -0.1  -0.1   0.0   0.0   0.1   516    1
## wyear[23]   0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.1  1200    1
## wyear[24]   0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.1  1200    1
## S           0.1     0.0 0.0   0.1   0.1   0.1   0.1   0.1   557    1
## U[1]        0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.2  1200    1
## U[2]        0.1     0.0 0.1   0.0   0.0   0.1   0.1   0.2  1200    1
## U[3]        0.0     0.0 0.1  -0.1   0.0   0.0   0.0   0.1  1200    1
## U[4]        0.0     0.0 0.1  -0.1  -0.1   0.0   0.0   0.1  1200    1
## U[5]        0.1     0.0 0.1   0.0   0.1   0.1   0.1   0.2  1200    1
## U[6]        0.1     0.0 0.1   0.0   0.0   0.1   0.1   0.2  1200    1
## U[7]        0.1     0.0 0.1  -0.1   0.0   0.1   0.1   0.2  1200    1
## U[8]        0.1     0.0 0.1   0.0   0.1   0.1   0.1   0.2  1200    1
## U[9]        0.1     0.0 0.1   0.0   0.0   0.1   0.1   0.2  1200    1
## U[10]       0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.1  1200    1
## U[11]       0.1     0.0 0.1   0.0   0.0   0.1   0.1   0.2  1200    1
## U[12]       0.0     0.0 0.1  -0.1   0.0   0.0   0.0   0.1  1200    1
## U[13]       0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.1  1200    1
## U[14]       0.0     0.0 0.1  -0.1  -0.1   0.0   0.0   0.1  1200    1
## U[15]       0.1     0.0 0.1   0.0   0.0   0.1   0.1   0.2  1200    1
## U[16]       0.1     0.0 0.1  -0.1   0.0   0.1   0.1   0.2  1200    1
## U[17]       0.1     0.0 0.1  -0.1   0.0   0.1   0.1   0.2  1200    1
## U[18]       0.0     0.0 0.1  -0.1   0.0   0.0   0.0   0.1  1200    1
## U[19]       0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.1  1200    1
## U[20]       0.0     0.0 0.1  -0.2  -0.1   0.0   0.0   0.1  1200    1
## U[21]       0.1     0.0 0.1   0.0   0.1   0.1   0.2   0.2  1200    1
## U[22]       0.2     0.0 0.1   0.0   0.1   0.2   0.2   0.3  1200    1
## U[23]       0.1     0.0 0.1   0.0   0.1   0.1   0.2   0.2  1200    1
## U[24]       0.2     0.0 0.1   0.1   0.1   0.2   0.2   0.3  1200    1
## U[25]       0.2     0.0 0.1   0.1   0.1   0.2   0.2   0.3  1200    1
## U[26]       0.1     0.0 0.1   0.0   0.0   0.1   0.1   0.2  1200    1
## U[27]       0.1     0.0 0.1   0.0   0.0   0.1   0.1   0.2  1200    1
## U[28]       0.1     0.0 0.1   0.0   0.0   0.1   0.1   0.2  1200    1
## U[29]       0.2     0.0 0.1   0.1   0.1   0.2   0.2   0.3  1200    1
## U[30]       0.1     0.0 0.1   0.0   0.1   0.1   0.2   0.3  1200    1
## U[31]       0.1     0.0 0.1   0.0   0.1   0.1   0.2   0.2  1200    1
## U[32]       0.1     0.0 0.1   0.0   0.1   0.1   0.2   0.2  1200    1
## U[33]       0.1     0.0 0.1   0.0   0.1   0.1   0.1   0.2  1200    1
## U[34]       0.1     0.0 0.1   0.0   0.1   0.1   0.2   0.2  1200    1
## U[35]       0.1     0.0 0.1   0.0   0.1   0.1   0.1   0.2  1200    1
## U[36]       0.1     0.0 0.1   0.0   0.1   0.1   0.2   0.2  1200    1
## U[37]       0.1     0.0 0.1  -0.1   0.0   0.1   0.1   0.2  1200    1
## U[38]       0.1     0.0 0.1   0.0   0.0   0.1   0.1   0.2  1200    1
## U[39]       0.1     0.0 0.1   0.0   0.0   0.1   0.1   0.2  1200    1
## U[40]       0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.2  1200    1
## U[41]       0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.1  1200    1
## U[42]      -0.1     0.0 0.1  -0.2  -0.1  -0.1   0.0   0.1  1200    1
## U[43]      -0.1     0.0 0.1  -0.2  -0.1  -0.1   0.0   0.1  1200    1
## U[44]      -0.1     0.0 0.1  -0.2  -0.1  -0.1  -0.1   0.0  1200    1
## U[45]       0.0     0.0 0.1  -0.1  -0.1   0.0   0.0   0.1  1200    1
## U[46]      -0.1     0.0 0.1  -0.2  -0.2  -0.1  -0.1   0.0  1200    1
## U[47]      -0.1     0.0 0.1  -0.2  -0.1  -0.1  -0.1   0.0  1200    1
## U[48]       0.0     0.0 0.1  -0.1  -0.1   0.0   0.0   0.1  1200    1
## U[49]      -0.1     0.0 0.1  -0.2  -0.1  -0.1  -0.1   0.0  1200    1
## U[50]      -0.2     0.0 0.1  -0.3  -0.2  -0.2  -0.1  -0.1  1200    1
## U[51]      -0.2     0.0 0.1  -0.3  -0.2  -0.2  -0.1   0.0  1200    1
## U[52]      -0.1     0.0 0.1  -0.2  -0.1  -0.1  -0.1   0.0  1200    1
## U[53]       0.1     0.0 0.1   0.0   0.0   0.1   0.1   0.2  1200    1
## U[54]       0.2     0.0 0.1   0.1   0.2   0.2   0.3   0.4  1200    1
## U[55]       0.1     0.0 0.1   0.0   0.1   0.1   0.2   0.3  1200    1
## U[56]       0.1     0.0 0.1   0.0   0.1   0.1   0.2   0.2  1200    1
## U[57]       0.1     0.0 0.1  -0.1   0.0   0.1   0.1   0.2  1200    1
## U[58]       0.0     0.0 0.1  -0.1   0.0   0.0   0.0   0.1  1200    1
## U[59]       0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.1  1200    1
## U[60]       0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.1  1200    1
## U[61]       0.0     0.0 0.1  -0.1   0.0   0.0   0.0   0.1  1200    1
## U[62]       0.1     0.0 0.1   0.0   0.0   0.1   0.1   0.2  1200    1
## U[63]       0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.2  1052    1
## U[64]       0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.1  1200    1
## U[65]       0.1     0.0 0.1   0.0   0.1   0.1   0.1   0.2  1200    1
## U[66]       0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.2  1200    1
## U[67]       0.0     0.0 0.1  -0.1   0.0   0.0   0.0   0.1  1200    1
## U[68]       0.1     0.0 0.1   0.0   0.0   0.1   0.1   0.2  1200    1
## U[69]       0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.1  1200    1
## U[70]       0.0     0.0 0.1  -0.1  -0.1   0.0   0.0   0.1  1200    1
## U[71]       0.0     0.0 0.1  -0.1  -0.1   0.0   0.0   0.1  1200    1
## U[72]      -0.1     0.0 0.1  -0.2  -0.1  -0.1  -0.1   0.0  1200    1
## U[73]      -0.1     0.0 0.1  -0.2  -0.1  -0.1   0.0   0.1   697    1
## U[74]      -0.1     0.0 0.1  -0.3  -0.2  -0.1  -0.1   0.0   865    1
## U[75]      -0.5     0.0 0.1  -0.6  -0.5  -0.5  -0.4  -0.3   567    1
## U[76]      -0.3     0.0 0.1  -0.5  -0.4  -0.3  -0.3  -0.2  1200    1
## U[77]       0.2     0.0 0.1   0.1   0.1   0.2   0.2   0.3  1200    1
## U[78]       0.2     0.0 0.1   0.1   0.1   0.2   0.2   0.3  1200    1
## U[79]       0.2     0.0 0.1   0.1   0.1   0.2   0.2   0.3  1200    1
## U[80]       0.1     0.0 0.1   0.0   0.1   0.1   0.2   0.3  1200    1
## U[81]       0.1     0.0 0.1  -0.1   0.0   0.1   0.1   0.2  1200    1
## U[82]       0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.1  1200    1
## U[83]       0.1     0.0 0.1   0.0   0.1   0.1   0.2   0.2  1200    1
## U[84]       0.1     0.0 0.1   0.0   0.1   0.1   0.2   0.3  1200    1
## U[85]       0.1     0.0 0.1  -0.1   0.0   0.1   0.1   0.2  1200    1
## U[86]      -0.1     0.0 0.1  -0.2  -0.1  -0.1   0.0   0.1  1200    1
## U[87]       0.0     0.0 0.1  -0.1   0.0   0.0   0.0   0.1  1200    1
## U[88]       0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.2  1200    1
## U[89]       0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.1  1200    1
## U[90]       0.0     0.0 0.1  -0.1   0.0   0.0   0.1   0.2  1200    1
## U[91]       0.1     0.0 0.1   0.0   0.0   0.1   0.1   0.2  1200    1
## U[92]       0.1     0.0 0.1   0.0   0.1   0.1   0.1   0.2  1200    1
## U[93]       0.1     0.0 0.1  -0.1   0.0   0.1   0.1   0.2  1200    1
## U[94]       0.2     0.0 0.1   0.0   0.1   0.2   0.2   0.3  1200    1
## U[95]       0.1     0.0 0.1  -0.1   0.0   0.1   0.1   0.2  1200    1
## lp__      153.9     0.3 4.6 143.7 151.0 154.3 157.4 161.5   288    1
## 
## Samples were drawn using NUTS(diag_e) at Thu Apr 24 22:58:01 2014.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).