#####################################################################################
#
# Application 3: Intradaily realized volatility of Alcoa stock (5m, 10m, 20m)
#
# Lopes and Tsay (2011)
# Particle Filters and Bayesian Inference in Financial Econometrics.
# Journal of Forecasting.
# 
#####################################################################################
#
# We entertain two realized volatility models: 
# (i) Three RV time series are modeled by independent univariate local level models
# (ii) Trivariate vector of RV time series is modeled by a multivariate local level
#
# Data: Intradaily realized volatility of Alcoa stock (5m, 10m, 20m) from 
# 2 January 2003 to 7 May 2004 for 340 observations.  The daily realized volatilities 
# used are the the sums of squares of intraday 5 min, 10 min and 20 min log returns 
# measured in percentages.
# Tsay (2005), Chapter 11: State-Space Models and Kalman Filter
# http://faculty.chicagobooth.edu/ruey.tsay/teaching/fts2/aa-3rv.txt
#
######################################################################################
#
# Author : Hedibert Freitas Lopes
#          Associate Professor of Econometrics and Statistics
#          The University of Chicago Booth School of Business
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoBooth.edu
#          URL: http://faculty.chicagobooth.edu/hedibert.lopes
#
#####################################################################################
rm(list=ls())
ldig=function(x,a,b){dgamma(1/x,a,b,log=TRUE)-2*log(x)}
skewness=function(x){n=length(x);s=sqrt(var(x));m=mean(x);return(sum(((x-m)/s)^3)*n/((n-1)*(n-2)))}
kurtosis=function(x){n=length(x);s=sqrt(var(x));m=mean(x);A=n*(n+1)/((n-1)*(n-2)*(n-3));B=3*(n-1)^2/((n-2)*(n-3));return(A*sum(((x-m)/s)^4)-B)}

PL.udlm = function(y,m0,C0,a0,b0,c0,d0,N){
  n     = length(y)
  mt    = rep(m0,N)
  Ct    = rep(C0,N)
  sig2  = 1/rgamma(N,a0,b0)
  tau2  = 1/rgamma(N,c0,d0)
  S     = matrix(c(a0,b0,c0,d0),N,4,byrow=T)
  sig2s = matrix(0,n,N)
  tau2s = matrix(0,n,N)
  xs    = matrix(0,n,N)
  like  = rep(0,n)
  for (t in 1:n){
  	  Wt   = tau2+sig2
  	  Rt   = Ct+tau2
  	  Qt   = Rt+sig2
  	  At   = Rt/Qt
  	  # Resampling s[t],sx[t],sig2,tau2
    w    = dnorm(y[t],mt,sqrt(Qt))
    like[t] = mean(w)
    k    = sample(1:N,size=N,prob=w,replace=T)
    mt   = mt[k];Ct=Ct[k]
    At   = At[k];Wt=Wt[k]
    Rt   = Rt[k];Qt=Qt[k]
    sig2 = sig2[k];tau2=tau2[k]
    S    = S[k,]
    # Sampling (x[t-1],x[t]) given (y^[t],tau2,sig2)
    V  = 1/(1/Ct+1/Wt)
    g  = V*(mt/Ct+y[t]/Wt)
    x1 = rnorm(N,g,sqrt(V))
    Cb = 1/(1/tau2+1/sig2)
    h  = Cb*(x1/tau2+y[t]/sig2)
    x  = rnorm(N,h,sqrt(Cb))
    # Updating state sufficient statistics
    mt = (1-At)*mt+At*y[t]
    Ct = Rt-At^2*Qt
    # Updating parameter sufficient statistics
    S[,1] = S[,1] + 0.5
    S[,2] = S[,2] + 0.5*(y[t]-x)^2
    S[,3] = S[,3] + 0.5
    S[,4] = S[,4] + 0.5*(x-x1)^2
    # Sampling parameters sig2 and tau2
    sig2  = 1/rgamma(N,S[,1],S[,2])
    tau2  = 1/rgamma(N,S[,3],S[,4])
    # Storing the output
    sig2s[t,] = sig2
    tau2s[t,] = tau2
    xs[t,]    = x
  }
  return(list(sig2s=sig2s,tau2s=tau2s,xs=xs,like=like))
}

PL.mdlm = function(y,c0,d0,eta0,S0,m0,C0,N){
  n      = nrow(y)
  p      = ncol(y)
  iS0    = solve(S0)
  x      = rnorm(N,m0,sqrt(C0))
  tau2   = 1/rgamma(N,c0,d0)
  Sigma  = array(0,c(N,p,p))
  iSigma = array(0,c(N,p,p))
  St     = array(0,c(N,p,p))
  Qt     = array(0,c(N,p,p))
  At     = matrix(0,N,p)
  for (j in 1:N){
    Sigma[j,,]  = rwishart(eta0,iS0)$IW
    iSigma[j,,] = solve(Sigma[j,,])
    St[j,,]     = S0
  }
  D       = matrix(1,p,p)
  mt      = rep(0,N)
  Ct      = rep(10,N)
  iW      = array(0,c(N,p,p))
  oneiW   = matrix(0,N,p)
  oneiSig = matrix(0,N,p)
  ct      = rep(c0,N)
  dt      = rep(d0,N)
  etat    = rep(eta0,N)
  draws   = array(0,c(n,N,8))
  w       = rep(0,N)
  like    = rep(0,n)
  for (t in 1:n){
  	  print(t)
    Rt = Ct+tau2
    for (j in 1:N){
      Qt[j,,]     = Rt[j]+Sigma[j,,]
    	  At[j,]      = matrix(Rt[j]*apply(solve(Qt[j,,]),2,sum),1,p)
      iW[j,,]     = solve(tau2[j]+Sigma[j,,])
      oneiW[j,]   = apply(iW[j,,],2,sum)
      oneiSig[j,] = apply(iSigma[j,,],2,sum) 
      w[j]        = dmnorm(y[t,],rep(mt[j],p),Qt[j,,],log=TRUE)
    }  
    like[t] = mean(exp(w))
    w       = exp(w-max(w))
    w       = w/sum(w)
    k       = sample(1:N,size=N,prob=w,replace=T)
    Sigma   = Sigma[k,,]
    iSigma  = iSigma[k,,]
    oneiSig = oneiSig[k,]
    iW      = iW[k,,]
    oneiW   = oneiW[k,]
    tau2    = tau2[k]
    Rt      = Rt[k]
    mt      = mt[k]
    Ct      = Ct[k]
    ct      = ct[k]
    dt      = dt[k]
    etat    = etat[k]
    St      = St[k,,]
    Qt      = Qt[k,,]
    At      = At[k,]
    V       = 1/(1/Ct+apply(iW,1,sum))
    g       = as.numeric(V*(mt/Ct+oneiW%*%y[t,]))
    x1      = rnorm(N,g,sqrt(V))
    Cbar    = 1/(1/tau2+apply(iSigma,1,sum))
    h       = Cbar*(x1/tau2+oneiSig%*%y[t,])
    x       = rnorm(N,h,sqrt(Cbar))
    ct      = ct+0.5
    dt      = dt+0.5*(x-x1)^2
    tau2    = 1/rgamma(N,ct,dt)
    etat    = etat+1
    for (j in 1:N){
    	  mt[j]       = mt[j]+At[j,]%*%(y[t,]-mt[j])
    	  Ct[j]       = Rt[j]-At[j,]%*%Qt[j,,]%*%At[j,]
      cross       = matrix(y[t,]-x[j],p,1)
      St[j,,]     = St[j,,]+cross%*%t(cross)
      iSt         = solve(St[j,,])
      draw        = rwishart(etat[j],iSt)
      Sigma[j,,]  = draw$IW
      iSigma[j,,] = draw$W
    }
    draws[t,,1] = tau2
    draws[t,,2] = x
    draws[t,,3] = Sigma[,1,1]
    draws[t,,4] = Sigma[,2,1]
    draws[t,,5] = Sigma[,2,2]
    draws[t,,6] = Sigma[,3,1]
    draws[t,,7] = Sigma[,3,2]
    draws[t,,8] = Sigma[,3,3]
  }
  return(list(draws=draws,like=like))
}

# Dataset
# -------
y1=c(3.3635,4.3848,1.5408,1.7437,107.0411,3.3061,5.7399,3.0006,2.9741,4.404,2.9325,2.5291,3.3107,4.1984,3.1586,6.8693,5.4874,
3.5684,7.1218,3.8051,3.6319,1.9308,4.9121,14.0632,8.8533,5.64,5.0164,4.9164,5.943,3.5707,5.7151,5.6497,3.1382,2.9794,7.9499,
3.2621,10.9732,5.0506,4.2809,4.4267,3.7911,1.6694,3.4518,3.0698,4.7244,3.3971,4.4795,5.4469,10.1704,8.1916,13.3029,4.0956,
5.0765,5.6913,6.7426,6.2963,5.189,2.8877,3.6595,3.3048,5.2429,6.2613,4.5942,4.2583,3.4015,4.8796,2.3002,3.7668,2.8985,3.723,
4.0437,2.4832,2.0142,1.8894,2.6252,1.7712,1.7464,2.8609,3.3228,1.7253,3.0478,3.7627,2.2816,2.9995,2.6456,2.5294,3.448,2.5817,
1.0716,1.8023,2.4987,1.7762,1.8979,1.2781,2.2896,3.6982,3.5916,1.4303,0.9581,2.2303,1.5319,5.3593,3.031,4.8053,3.1961,2.6163,
3.3775,5.7123,2.7711,3.7482,2.1212,1.9637,1.9014,1.5668,4.7035,2.5363,3.2273,1.8119,1.2525,2.2765,1.7602,1.7115,1.8768,3.1587,
2.0543,1.6741,6.3634,1.7686,1.7783,5.8654,2.0188,1.8085,1.6768,3.2024,3.5515,3.8947,1.9297,3.4713,2.4262,2.725,2.764,2.7329,
3.8417,3.1981,2.1192,3.319,4.5216,3.6186,2.6966,3.3665,2.1307,1.3193,1.7938,1.7457,1.4514,1.8184,1.3031,1.6634,2.2056,1.9387,
2.937,1.9929,4.3651,2.5257,1.0882,2.3673,2.2522,2.8816,1.7091,2.1168,1.4164,2.0983,3.2058,4.0238,1.6034,5.1311,2.2936,1.8926,
2.143,1.4239,2.6464,1.6464,2.8841,3.0977,4.6001,3.3997,1.9265,5.2172,1.6076,2.8704,6.2519,1.697,1.8732,2.3129,1.9364,1.5129,
0.7924,2.1521,1.709,2.9699,1.0849,1.0892,0.9484,1.8906,2.3129,1.0992,1.6592,2.9976,2.1879,2.2682,2.9851,3.3353,1.6794,2.981,
1.5941,2.3089,2.5672,1.8663,2.754,1.3907,3.6308,4.087,3.561,1.6556,2.8988,1.379,4.8036,1.8944,2.8527,0.4795,1.9401,0.8485,
2.0725,1.5146,2.3422,0.964,2.4085,2.3175,1.9954,0.7189,1.4913,1.3322,1.6372,2.1888,4.6134,2.942,2.1492,1.4034,1.9647,1.6742,
1.8931,1.0535,2.5975,1.4315,2.6529,1.9291,3.325,6.4407,3.3292,3.6127,2.696,3.2308,2.7415,1.9109,2.7028,2.0763,3.6868,1.4607,
1.2173,2.8695,3.8392,4.1788,2.7565,2.4172,2.5231,3.3259,1.98,1.283,1.5097,2.4873,3.2211,2.083,1.2458,1.0851,2.1153,1.8359,
3.3339,2.1944,1.3395,2.7383,4.6535,2.6592,2.389,2.8173,1.9898,3.1297,1.4581,2.272,2.7152,4.3215,2.0795,2.8156,2.2521,2.1545,
1.7163,2.3143,2.8758,2.7332,3.8024,2.0057,1.3821,1.2844,0.9123,1.5872,2.5897,1.7687,1.8486,1.748,13.3502,3.109,1.3721,3.107,
2.4874,3.2602,1.6436,1.7224,1.8612,17.7923,5.7749,3.768,5.1514,4.7074,5.5567,10.8065,3.2743,3.0118,4.8383,1.4688,3.5661,3.6526)

y2=c(3.4745,4.147,1.2124,1.8723,11.0797,3.2015,5.8178,2.194,2.597,5.926,3.3583,3.1004,3.1949,3.0953,3.0546,5.8683,4.9194,
2.9026,5.9915,3.6607,2.7014,1.8549,5.7454,3.8358,9.3412,4.7918,5.0295,5.1436,6.1604,3.2566,4.7144,4.4854,2.5977,4.0665,
8.8823,3.4285,5.2562,4.1971,2.7653,4.4831,3.7528,2.1363,3.0092,3.4804,2.5672,3.1837,4.5751,4.5726,3.4218,8.3076,15.573,
3.1188,3.3909,4.1333,6.7761,6.6053,3.9847,2.3069,2.3601,3.5887,6.1575,8.7782,3.0794,3.4978,3.2078,3.9,2.5396,3.7344,2.401,
4.4359,2.0062,1.8092,2.0674,1.5805,2.5346,1.4713,1.6664,3.7704,3.1881,1.4827,2.6197,4.5255,2.3745,2.2796,1.9674,1.9773,
2.6925,2.6341,0.8949,1.3988,3.1687,1.7112,2.0132,1.3439,1.9521,3.9489,3.619,1.8653,1.0237,2.4006,1.5105,6.1684,2.7774,
4.5807,2.7089,3.1904,2.7875,6.4343,3.2695,2.1194,2.1098,1.2263,2.3053,1.4845,4.5845,2.503,2.3832,1.9295,1.5617,1.421,1.9033,
1.7105,2.2195,2.7157,2.1963,1.543,5.3884,1.8786,1.9635,3.8441,1.8782,1.5463,1.5694,3.3205,1.7537,3.0671,1.6297,4.3007,2.7463,
2.0095,3.8564,3.5055,3.9339,2.1294,2.2021,4.2113,5.002,2.8519,2.2789,3.3877,1.3113,1.3841,1.707,2.1297,1.4556,1.7629,1.1851,
1.7019,2.2483,1.6893,2.7568,1.969,3.4255,3.0274,1.1874,2.7625,2.5329,3.6644,1.4738,1.8363,1.1711,2.7982,4.273,1.8007,1.659,
7.7296,2.298,1.8201,1.7699,1.1194,3.8116,1.456,2.3592,3.2543,5.7432,3.2952,2.3301,5.8647,1.0407,2.927,9.5184,1.4462,1.9724,
2.3429,2.2228,1.739,0.8935,1.666,1.3037,2.4128,1.1326,0.8942,1.0162,2.1681,1.4579,0.8445,1.2125,3.1824,2.5537,1.7977,3.6401,
4.6919,1.93,2.0406,1.2317,2.4069,2.7515,2.3955,2.8676,1.3392,5.168,4.3562,3.3759,1.4837,3.5352,0.7982,5.8114,2.1643,1.6724,
0.3007,2.1877,0.8406,1.7749,1.696,2.021,0.7428,2.1738,2.2588,1.1623,0.5132,1.6509,1.4135,1.6797,2.5426,3.4689,2.3829,1.6848,
1.8694,1.111,1.9962,1.4167,1.2138,2.6778,1.1748,3.0136,2.1112,2.9758,8.0769,4.0148,3.9869,2.314,2.3194,3.4218,1.1506,1.7263,
2.1079,2.9261,1.2024,1.156,3.2218,4.4327,5.4168,3.6871,1.9243,1.9171,3.2805,2.27,0.918,1.7007,3.0829,3.352,1.8059,1.2941,
0.998,1.6861,2.0414,2.289,1.805,1.0193,3.1129,6.2269,2.2281,2.3667,1.2504,2.6165,1.855,1.7612,2.5231,3.2313,5.2266,2.1601,
2.8165,2.2392,1.2924,1.3721,1.7077,4.1317,1.9789,5.6926,1.7661,0.837,0.6544,0.7281,1.8337,1.406,1.6339,1.8427,1.358,1.7583,
4.441,1.7173,3.1552,2.5848,3.3807,1.8144,1.6676,1.3438,19.1777,5.4538,3.1398,3.3128,3.4229,8.3439,3.7491,3.0383,2.5793,7.1472,
1.4514,4.2928,3.5175)

y3 = c(3.0692,6.7121,0.5711,2.122,14.5181,3.303,2.765,1.9447,1.4474,6.5017,2.7507,3.7948,3.7756,3.1704,3.4981,6.8406,3.7551,
4.8222,3.9653,2.658,3.3734,1.0804,6.7529,1.7216,5.5931,2.7267,4.7233,4.2175,2.5537,1.6059,6.7449,6.5552,2.6176,2.0435,11.3686,
2.9554,3.7618,4.6676,2.3042,2.9324,2.8326,2.2052,3.0339,4.9699,1.7451,2.7317,5.5938,3.7813,3.2598,7.8057,12.5508,3.4027,2.3163,
4.9057,6.3096,6.3244,3.4566,2.4782,2.792,2.0617,6.1016,7.5055,3.3044,3.6158,2.5028,4.6102,3.0912,4.0922,2.6161,5.1447,1.852,
2.376,1.5587,1.7993,2.0181,1.3669,1.8802,4.2498,1.5221,1.2279,1.9778,3.429,2.7599,3.483,1.9109,1.5576,2.6136,2.3325,0.6791,
0.6751,3.3529,2.3148,2.128,1.9298,1.6712,2.7863,3.6446,2.4316,1.2834,3.1579,1.305,6.4936,1.8195,4.4878,3.0601,3.3746,2.8478,
8.1738,4.1288,1.5081,1.1789,1.2877,3.2212,1.4368,2.989,2.0541,2.632,1.4976,1.8638,1.1949,2.0815,1.3515,2.6168,2.5188,3.1049,
1.586,1.8715,2.388,1.6246,3.8855,1.8175,1.9876,1.4315,3.8023,1.6497,1.8234,1.2766,3.2726,3.3855,2.3532,5.1794,4.9593,2.9302,
2.6604,2.0022,5.2051,3.4881,2.1066,2.3471,2.1273,1.3812,1.0161,1.3418,2.5714,1.6548,2.4974,0.6987,1.9332,2.1974,0.953,2.6124,
2.2209,4.6254,1.9715,1.814,2.9249,2.6555,4.6945,1.0214,0.925,1.2941,4.2612,3.7743,2.7033,1.1286,4.3708,2.7175,2.6416,1.22,
1.9848,4.4325,1.5964,3.4269,4.8462,5.2844,3.0862,3.3093,6.9519,1.4129,3.5558,6.4416,1.9195,2.0843,2.0544,2.7649,1.4205,0.7672,
1.3993,1.3193,2.172,1.7648,0.4611,1.2857,2.3191,1.0733,0.7451,0.5596,4.4416,2.8547,1.527,3.107,3.8656,2.5505,2.5417,1.3014,
1.9364,3.3863,2.6213,3.4443,1.6597,6.1557,2.9505,1.7269,1.3044,4.7293,0.6855,6.5098,1.8704,2.0294,0.4252,2.6016,0.7381,2.3344,
2.2122,2.0229,0.5458,2.2005,2.0219,0.9058,0.2669,2.0008,1.4464,0.9676,2.0688,2.3751,3.9783,1.5709,1.3658,6.2888,1.8073,1.1273,
1.2406,1.7327,1.3453,3.809,1.1641,1.9583,9.6052,4.0337,4.3449,1.4234,2.5998,4.7903,1.4189,2.307,2.0939,2.264,1.7613,1.0687,
5.5588,4.9795,3.2766,5.3183,1.1397,1.3891,3.3294,1.0457,0.8432,0.9447,1.9117,2.3126,1.5158,1.4965,1.0103,1.3914,2.3527,1.629,
2.1664,0.6762,3.0233,4.6574,2.8416,2.1978,1.1459,1.6478,1.489,1.2585,2.4157,3.1367,3.7711,1.8246,2.4722,1.861,1.3226,1.2434,
1.6485,4.3189,1.3671,5.5476,1.3963,0.4982,0.6758,0.4957,1.5065,1.4978,1.3344,1.2167,1.2679,1.3033,3.8509,1.05,3.1512,2.6406,
3.3868,1.5142,2.069,1.6506,28.1534,5.6816,3.4632,3.073,1.8359,9.9087,3.7128,1.7533,3.0946,2.5901,1.1785,3.0513,4.8249)

y = cbind(y1,y2,y3)
y     = log(data)
n     = nrow(y) 
p     = ncol(y)
names = c("5-minute","10-minute","20-minute")

pdf(file="rv-timeseries.pdf",width=8,height=12)
par(mfrow=c(3,2))
i=1;ts.plot(y[,i],xlab="Day",ylab="",main=names[i],ylim=range(y))
plot(y[,1:2],xlab=names[1],ylab=names[2],pch=16,xlim=range(y),ylim=range(y));abline(0,1)
i=2;ts.plot(y[,i],xlab="Day",ylab="",main=names[i],ylim=range(y))
plot(y[,c(1,3)],xlab=names[1],ylab=names[3],pch=16,xlim=range(y),ylim=range(y));abline(0,1)
i=3;ts.plot(y[,i],xlab="Day",ylab="",main=names[i],ylim=range(y))
plot(y[,2:3],xlab=names[2],ylab=names[3],pch=16,xlim=range(y),ylim=range(y));abline(0,1)
dev.off()

round(cbind(apply(y,2,mean),apply(y,2,median),apply(y,2,skewness),apply(y,2,kurtosis)),3)
#5-min   0.992  0.977 1.091 5.479
#10-min  0.913  0.871 0.153 0.769
#20-min  0.850  0.847 0.034 0.843

round(var(y),3)
#[1,] 0.314 0.270 0.258
#[2,] 0.270 0.317 0.307
#[3,] 0.258 0.307 0.396

round(cor(y),3)
#[1,] 1.000 0.857 0.732
#[2,] 0.857 1.000 0.865
#[3,] 0.732 0.865 1.000

# Particle learning (Univariate)
# ------------------------------
set.seed(13579)
N     = 10000
m0    = 0.0
C0    = 1.0
a0    = 10
b0    = (a0+1)*0.10
c0    = 10
d0    = (c0+1)*0.05
c(b0/(a0-1),b0/(a0+1),1/qgamma(0.975,a0,b0),1/qgamma(0.025,a0,b0))
c(d0/(c0-1),d0/(c0+1),1/qgamma(0.975,c0,d0),1/qgamma(0.025,c0,d0))
qu    = array(0,c(3,3,n,3))
sig2u = array(0,c(3,n,N))
tau2u = array(0,c(3,n,N))
xu    = array(0,c(3,n,N))
like  = matrix(0,3,n)
for (i in 1:3){
	print(i)
  plu        = PL.udlm(y[,i],m0,C0,a0,b0,c0,d0,N)
  sig2u[i,,] = plu$sig2s
  tau2u[i,,] = plu$tau2s
  xu[i,,]    = plu$xs
  like[i,]   = plu$like
  qu[i,1,,]  = t(apply(sig2u[i,,],1,quantile,c(0.025,0.5,0.975)))
  qu[i,2,,]  = t(apply(tau2u[i,,],1,quantile,c(0.025,0.5,0.975)))
  qu[i,3,,]  = t(apply(xu[i,,],1,quantile,c(0.025,0.5,0.975)))
}

pdf(file="rv-univariate1.pdf",width=8,height=12)
par(mfrow=c(3,2))
ind = 50:n
for (i in 1:3){
  ts.plot(qu[i,1,,],xlab="Time",ylab=expression(sigma^2),ylim=range(qu[,1,ind,]))
  title(names[i])
  ts.plot(qu[i,2,,],xlab="Time",ylab=expression(tau^2),ylim=range(qu[,2,ind,]))
  title(names[i])
}
dev.off()

pdf(file="rv-univariate2.pdf",width=8,height=12)
par(mfrow=c(3,2))
U       = c(25,75)
breaks1 = seq(min(sig2u[,n,]),max(sig2u[,n,]),length=30)
breaks2 = seq(min(tau2u[,n,]),max(tau2u[,n,]),length=30)
for (i in 1:3){
  hist(sig2u[i,n,],xlim=range(sig2u[,n,]),xlab=expression(sigma^2),
       ylab="",prob=TRUE,main=names[i],ylim=c(0,U[1]),breaks=breaks1)
  box()
  hist(tau2u[i,n,],xlim=range(tau2u[,n,]),xlab=expression(tau^2),
       ylab="",prob=TRUE,main=names[i],ylim=c(0,U[2]),breaks=breaks2)
  box()
}
dev.off()

pdf(file="rv-univariate3.pdf",width=8,height=12)
par(mfrow=c(3,1))
cols = c(grey(0.6),1,grey(0.6))
for (i in 1:3)
  ts.plot(qu[i,3,ind,],ylim=range(qu[,3,ind,]),main=names[i],ylab="",col=cols,lwd=3)
dev.off()

# Particle learning (Multivariate)
# --------------------------------
set.seed(13579)
N    = 2000
eta0 = 2*a0
S0   = diag(2*b0,p)
S0/(eta0-p-1)
S0/(eta0+p+1)
plm  = PL.mdlm(y,c0,d0,eta0,S0,m0,C0,N)
likem = plm$like
plm  = plm$draws
qm   = array(0,c(n,8,3))
qc   = array(0,c(n,3,3))
cs   = array(0,c(n,N,3))
cs[,,1] = plm[,,4]/sqrt(plm[,,3]*plm[,,5])
cs[,,2] = plm[,,6]/sqrt(plm[,,3]*plm[,,8])
cs[,,3] = plm[,,7]/sqrt(plm[,,5]*plm[,,8])
for (i in 1:8){
  qm[,i,1] = apply(plm[,,i],1,quantile,0.025)
  qm[,i,2] = apply(plm[,,i],1,quantile,0.5)
  qm[,i,3] = apply(plm[,,i],1,quantile,0.975)
}
for (i in 1:3){
  qc[,i,1] = apply(cs[,,i],1,quantile,0.025)
  qc[,i,2] = apply(cs[,,i],1,quantile,0.5)
  qc[,i,3] = apply(cs[,,i],1,quantile,0.975)
}  

likeu = apply(like,2,prod)
llikeu = log(likeu)
llikem = log(likem)

BF = likem/likeu
par(mfrow=c(1,1))
plot(BF,ylim=c(0,20),type="h",xlab="Time",ylab="Predictive ratio")
abline(h=1,lwd=2)

par(mfrow=c(1,1))
plot(log(cumprod(BF)),type="h",xlab="Time",ylab="Log cumulative Bayes factor")
abline(h=1,lwd=2)

# Output
# ------
pdf(file="rv-multivariate1.pdf",width=8,height=8)
ind = 50:n
par(mfrow=c(3,3))
cols = c(grey(0.5),1,grey(0.5))
ts.plot(qm[,3,],col=cols,ylim=range(qm[ind,3,]),xlab="",ylab="",main=expression(sigma[1]^2))
ts.plot(qc[,1,],col=cols,ylim=range(qc[ind,1,]),xlab="",ylab="",main=expression(rho[12]))
ts.plot(qc[,2,],col=cols,ylim=range(qc[ind,2,]),xlab="",ylab="",main=expression(rho[13]))
ts.plot(qm[,4,],col=cols,ylim=range(qm[ind,4,]),xlab="",ylab="",main=expression(sigma[21]))
ts.plot(qm[,5,],col=cols,ylim=range(qm[ind,5,]),xlab="",ylab="",main=expression(sigma[2]^2))
ts.plot(qc[,3,],col=cols,ylim=range(qc[ind,3,]),xlab="",ylab="",main=expression(rho[23]))
ts.plot(qm[,6,],col=cols,ylim=range(qm[ind,6,]),xlab="",ylab="",main=expression(sigma[31]))
ts.plot(qm[,7,],col=cols,ylim=range(qm[ind,7,]),xlab="",ylab="",main=expression(sigma[32]))
ts.plot(qm[,8,],col=cols,ylim=range(qm[ind,8,]),xlab="",ylab="",main=expression(sigma[3]^2))
dev.off()

pdf(file="rv-multivariate2.pdf",width=8,height=12)
par(mfrow=c(3,1))
ts.plot(qm[,1,],ylim=range(qm[ind,1,]),ylab=expression(tau^2),main="(a)")
hist(plm[n,,1],breaks=seq(min(plm[n,,1]),max(plm[n,,1]),length=30),xlab="",ylab="",prob=TRUE,main="(b)")
ts.plot(qm[,2,],ylim=range(qu[,3,ind,],qm[ind,2,]),xlab="",lwd=3,ylab=expression(x[t]),col=cols,main="(c)")
dev.off()

pdf(file="rv-multivariate3.pdf",width=8,height=12)
par(mfrow=c(2,1))
plot(density(tau2u[1,n,]),xlab="",ylab="",xlim=c(0.015,0.075),
     main=expression(tau^2),lwd=2,col=grey(0.5),ylim=c(0,100),lty=1)
lines(density(tau2u[2,n,]),lwd=2,col=grey(0.5),lty=2)
lines(density(tau2u[3,n,]),lwd=2,col=grey(0.5),lty=3)
lines(density(plm[n,,1]),lwd=3)
plot(density(sig2u[1,n,]),xlab="",ylab="",col=grey(0.5),xlim=c(0.1,0.4),
     main=expression(sigma^2),lwd=2,ylim=c(0,30))
lines(density(plm[n,,3]),lwd=3)
lines(density(sig2u[2,n,]),lwd=2,col=grey(0.5),lty=2)
lines(density(plm[n,,5]),lwd=3,lty=2)
lines(density(sig2u[3,n,]),lwd=2,col=grey(0.5),lty=3)
lines(density(plm[n,,8]),lwd=3,lty=3)
text(0.13,25,"5m",cex=1.5)
text(0.21,21,"10m",cex=1.5)
text(0.31,16,"20m",cex=1.5)
dev.off()

pdf(file="rv-multivariate4.pdf",width=8,height=12)
par(mfrow=c(3,2))
breaks = seq(0,1,length=50)
ts.plot(qc[,1,],col=cols,ylim=c(0,1),xlab="Time",ylab=expression(rho[12]),main="")
hist(cs[n,,1],prob=TRUE,breaks=breaks,xlab=expression(rho[12]),ylab="",main="",ylim=c(0,17))
ts.plot(qc[,2,],col=cols,ylim=c(0,1),xlab="Time",ylab=expression(rho[13]),main="")
hist(cs[n,,2],prob=TRUE,breaks=breaks,xlab=expression(rho[13]),ylab="",main="",ylim=c(0,17))
ts.plot(qc[,3,],col=cols,ylim=c(0,1),xlab="Time",ylab=expression(rho[23]),main="")
hist(cs[n,,3],prob=TRUE,breaks=breaks,xlab=expression(rho[23]),ylab="",main="",ylim=c(0,17))
dev.off()

pdf(file="rv-multivariate5.pdf",width=12,height=14)
par(mfrow=c(2,1))
ini = c(1,171)
fin = c(170,n)
half = c("Time (first half)","Time (second half)")
for (i in 1:2){
  ind = ini[i]:fin[i]
  ind1 = trunc(seq(ind[1],ind[170],length=5))
  plot(ind,qu[1,3,ind,2],ylim=range(qu[,3,ind,2],qm[ind,2,2]),xlab=half[i],main="",ylab=expression(x[t]),type="l",axes=FALSE,lwd=2)
  axis(2);box();axis(1,at=ind1,lab=ind1)
  lines(ind,qu[2,3,ind,2],,lwd=2,lty=2)
  lines(ind,qu[3,3,ind,2],lwd=2,lty=3)
  lines(ind,qm[ind,2,2],lwd=4)
  if (i==1){
    legend(100,2.45,legend=c(names,"multivariate"),lty=c(1,2,3,1),lwd=c(2,2,2,4),bty="n",cex=2)
  }
}
dev.off()