| | 140 | // Initialize volatility shocks by correlating the current (supplied) set of |
|---|
| | 141 | // iid random-walk variates |
|---|
| | 142 | |
|---|
| | 143 | // for each multivariate sample... |
|---|
| | 144 | for(k1=0; k1<Nz[1]; k1++) { |
|---|
| | 145 | |
|---|
| | 146 | // For each time series... |
|---|
| | 147 | for(k2=0; k2<Nz[0]; k2++) { |
|---|
| | 148 | // The matrix product of the cross-correlations and the independent |
|---|
| | 149 | // variates for this sample... |
|---|
| | 150 | sum = 0; |
|---|
| | 151 | for(km=0; km<Nz[0]; km++) |
|---|
| | 152 | sum += RV2(k2, km) * Y2(k2, 2); |
|---|
| | 153 | // ...is the correlated volatility shock for this time series and sample |
|---|
| | 154 | V2(k2, k1) = sum; |
|---|
| | 155 | } |
|---|
| | 156 | } |
|---|
| | 157 | |
|---|
| | 158 | |
|---|
| | 159 | // Hybrid Gibbs sampler with Metropolis-Hastings accept/reject step for each |
|---|
| | 160 | // conditional density draw. |
|---|
| | 161 | |
|---|
| 165 | | // Proceed with the random walk only if this isn't the first, |
|---|
| 166 | | // initialization iteration... |
|---|
| 167 | | if(k0 > 0) { |
|---|
| 168 | | // For each oversampling iteration... |
|---|
| 169 | | for(k3=0; k3<n3; k3++) { |
|---|
| 170 | | // Generate a uniform, symmetric, random walk proposal |
|---|
| 171 | | normp = Y2(k2, 2) + wiggle * (drand48() - 0.5); |
|---|
| 172 | | // Do the ratio in log space |
|---|
| 173 | | L1 = -0.5 * pow(normp, 2); |
|---|
| 174 | | dL = L1 - L0; |
|---|
| 175 | | // The uniform random variate thus needs to be in logspace as well, |
|---|
| 176 | | // but it's not always even needed. Thus the log operation is done |
|---|
| 177 | | // just once, and only some of the time. |
|---|
| 178 | | if(dL > 0 || log(drand48()) < dL) { |
|---|
| 179 | | // Update current value and its log probability when proposal |
|---|
| 180 | | // accepted |
|---|
| 181 | | Y2(k2, 2) = normp; |
|---|
| 182 | | L0 = L1; |
|---|
| 183 | | } |
|---|
| | 179 | // For each oversampling iteration... |
|---|
| | 180 | for(k3=0; k3<n3; k3++) { |
|---|
| | 181 | // Generate a uniform, symmetric, random walk proposal |
|---|
| | 182 | normp = Y2(k2, 2) + wiggle * (drand48() - 0.5); |
|---|
| | 183 | // Do the ratio in log space |
|---|
| | 184 | L1 = -0.5 * pow(normp, 2); |
|---|
| | 185 | dL = L1 - L0; |
|---|
| | 186 | // The uniform random variate thus needs to be in logspace as well, but |
|---|
| | 187 | // it's not always even needed. Thus the log operation is done just |
|---|
| | 188 | // once, and only some of the time. |
|---|
| | 189 | if(dL > 0 || log(drand48()) < dL) { |
|---|
| | 190 | // Update current value and its log probability when proposal |
|---|
| | 191 | // accepted |
|---|
| | 192 | Y2(k2, 2) = normp; |
|---|
| | 193 | L0 = L1; |
|---|
| 202 | | // Now do the VAR(1) to get a multivariate log-volatility proposal. The |
|---|
| 203 | | // proposals are generated from "left to right" based solely on the |
|---|
| 204 | | // random-walk normal proposals, so there's no need to worry about |
|---|
| 205 | | // overwriting anything. |
|---|
| 206 | | |
|---|
| 207 | | // For each time series... |
|---|
| 208 | | for(k2=0; k2<Nz[0]; k2++) { |
|---|
| 209 | | if(k1 == 0) { |
|---|
| 210 | | // The log-volatility value of the first multivariate sample is just the |
|---|
| 211 | | // offset scaled by its VAR gain. |
|---|
| 212 | | H2(k2, 0) = D1(k2) / (1 - E2(k2, k2)); |
|---|
| 213 | | } else { |
|---|
| 214 | | // For all other samples, the offset plus the current proposed shock... |
|---|
| 215 | | sum = D1(k2) + Y2(k2, 3); |
|---|
| 216 | | // ...plus the dot product for the VAR(1) term... |
|---|
| 217 | | for(km=0; km<Nz[0]; km++) |
|---|
| 218 | | sum += E2(k2, km) * H2(k2, k1-1); |
|---|
| 219 | | // ...is the proposed log-volatility value |
|---|
| 220 | | H2(k2, k1) = sum; |
|---|
| 221 | | } |
|---|
| | 212 | // Now do the VAR(1) starting with this proposed volatility shock and |
|---|
| | 213 | // proceeding through the other shocks subsequent to it to compute |
|---|
| | 214 | // Pr(hp[t]|x,h[:t],h[t+1:]). |
|---|
| | 215 | L = 0; |
|---|
| | 216 | if(k3+nv >= Nz[1]) { |
|---|
| | 217 | nd = Nz[1]; |
|---|
| | 218 | } else { |
|---|
| | 219 | nd = k3+nv; |
|---|
| 223 | | |
|---|
| 224 | | // Now inverse-scale the innovations for this sample by the square root of |
|---|
| 225 | | // the volatilities in preparation for decorrelating the innovations |
|---|
| 226 | | for(k2=0; k2<Nz[0]; k2++) |
|---|
| 227 | | Y2(k2, 4) = exp(-0.5 * H2(k2, k1)) * X2(k2, k1); |
|---|
| 228 | | |
|---|
| 229 | | // Now decorrelate the equi-variance innovations in preparation for the |
|---|
| 230 | | // log-likelihood computation for this multivariate log-volatility proposal |
|---|
| 231 | | |
|---|
| 232 | | // For each time series... |
|---|
| 233 | | for(k2=0; k2<Nz[0]; k2++) { |
|---|
| 234 | | // The product of the inverted cross-correlation matrix and the |
|---|
| 235 | | // correlated variates... |
|---|
| 236 | | sum = 0; |
|---|
| 237 | | for(km=0; km<Nz[0]; km++) |
|---|
| 238 | | sum += RII2(k2, km) * Y2(k2, 4); |
|---|
| 239 | | // ...is the decorrelated multivariate innovation for this sample, given |
|---|
| 240 | | // its proposed log-volatility |
|---|
| 241 | | Y2(k2, 5) = sum; |
|---|
| 242 | | } |
|---|
| 243 | | |
|---|
| 244 | | // Now compute log Pr(x[t]|hp[t]), the log-likelihood of the decorrelated |
|---|
| 245 | | // multivariate innovation, conditional on the proposed log-volatility |
|---|
| 246 | | // hp[t], which has been generated conditional on the log-volatility of the |
|---|
| 247 | | // previous sample. Consider the proposal as a draw from Pr(hp[t]|h[t-1]). |
|---|
| 248 | | |
|---|
| 249 | | sum = 0; |
|---|
| 250 | | // For each time series... |
|---|
| 251 | | for(k2=0; k2<Nz[0]; k2++) { |
|---|
| 252 | | // Offset the decorrelated innovation by the mean, which is the |
|---|
| 253 | | // correlation between the normal variates for innovation and proposed |
|---|
| 254 | | // volatility shocks, times the value of the normal variate value for the |
|---|
| 255 | | // proposed volatility shock for this sample |
|---|
| 256 | | Y2(k2, 5) = Y2(k2, 5) - G1(k2) * H2(k2, k1); |
|---|
| 257 | | // The log-likelihood is simply the argument of the exponential of the |
|---|
| 258 | | // normal pdf, with the reciprocal scaling of the exponential being done |
|---|
| 259 | | // in log space by subtraction. |
|---|
| 260 | | sum -= pow(Y2(k2, 5) / Y2(k2, 0), 2) + Y2(k2, 1); |
|---|
| 261 | | // Since the innovation has been inverse-scaled by the square root of the |
|---|
| 262 | | // proposed volatility, the probability density needs to be scaled as |
|---|
| 263 | | // well. In log space, this is a subtraction of half the proposed |
|---|
| 264 | | // log-volatility. |
|---|
| 265 | | sum -= 0.5 * H2(k2, k1); |
|---|
| 266 | | } |
|---|
| 267 | | P2(0, k1) = sum; |
|---|
| 268 | | |
|---|
| 269 | | // The probability density that we need for our Metropolis-Hastings |
|---|
| 270 | | // accept/reject decision is Pr(hp[t]|x[t],h[t+1]), where h[t+1] is the |
|---|
| 271 | | // current log-volatility of the next sample. |
|---|
| 272 | | // |
|---|
| 273 | | // According to Bayes' Rule, the required density is proportional to |
|---|
| 274 | | // Pr(x[t]|hp[t],h[t+1]) * Pr(hp[t]|h[t+1]). So, we need to compute and |
|---|
| 275 | | // scale by the probability density of the log-volatility proposal given |
|---|
| 276 | | // the current log-volatility of the next sample, which is, by Bayes rule: |
|---|
| 277 | | // |
|---|
| 278 | | // Pr(h[t+1]|hp[t]) * Pr(hp[t]) / Pr(h[t+1]) |
|---|
| 279 | | // <-----+--------------> |
|---|
| 280 | | // | |
|---|
| 281 | | // +--- But this = 1 because both are equiprobable |
|---|
| 282 | | // |
|---|
| 283 | | // So, we now compute Pr(h[t+1]|hp[t]), which is Pr(h[t+1] - E2*hp[t]) |
|---|
| 284 | | // unless we're at the last sample, in which case it is = 1. |
|---|
| 285 | | |
|---|
| 286 | | if(k1 < Nz[1]) { |
|---|
| 287 | | // First we need to remove the offset and the current proposal's VAR term |
|---|
| 288 | | // from the next sample's current log-volatility |
|---|
| 289 | | |
|---|
| | 221 | |
|---|
| | 222 | // For the current sample and some samples subsequent to it... |
|---|
| | 223 | for(k3=k1; k3<nd; k3++) { |
|---|
| | 224 | // We alternate between the two columns of the log-volatility scratchpad |
|---|
| | 225 | // for current & previous multivariate log-volatility values |
|---|
| | 226 | kh = k3 % 2; |
|---|
| 292 | | // ...the offset |
|---|
| 293 | | sum = D1(k2); |
|---|
| 294 | | // ...plus the matrix product of the Lag-1 cross-correlations and the |
|---|
| 295 | | // proposed multivariate log-volatility components... |
|---|
| 296 | | for(km=0; km<Nz[0]; km++) |
|---|
| 297 | | sum += E2(k2, km) * H2(k2, k1); |
|---|
| 298 | | // ...is the VAR term to subtract from the next sample's log-volatility |
|---|
| 299 | | // value |
|---|
| 300 | | Y2(k2, 6) = H2(k2, k1+1) - sum; |
|---|
| 301 | | } |
|---|
| 302 | | |
|---|
| 303 | | // Then we need to decorrelate the residuals to get independent normal |
|---|
| 304 | | // variates |
|---|
| 305 | | |
|---|
| | 229 | if(k1 == 0) { |
|---|
| | 230 | // The log-volatility value of the first multivariate sample is just |
|---|
| | 231 | // the offset scaled by its VAR gain, plus the proposed shock |
|---|
| | 232 | H2(k2, 0) = D1(k2) / (1 - E2(k2, k2)) + Y2(k2, 3); |
|---|
| | 233 | } else { |
|---|
| | 234 | // For all other samples, the offset plus the proposed shock... |
|---|
| | 235 | sum = D1(k2) + Y2(k2, 3); |
|---|
| | 236 | // ...plus the dot product for the VAR(1) term... |
|---|
| | 237 | for(km=0; km<Nz[0]; km++) |
|---|
| | 238 | sum += E2(k2, km) * H2(k2, 1-kh); |
|---|
| | 239 | // ...is the proposed log-volatility value |
|---|
| | 240 | H2(k2, kh) = sum; |
|---|
| | 241 | } |
|---|
| | 242 | } |
|---|
| | 243 | |
|---|
| | 244 | // Now inverse-scale the innovations for this sample by the square root |
|---|
| | 245 | // of the volatilities in preparation for decorrelating the innovations |
|---|
| | 246 | for(k2=0; k2<Nz[0]; k2++) |
|---|
| | 247 | Y2(k2, 4) = exp(-0.5 * H2(k2, kh)) * X2(k2, k3); |
|---|
| | 248 | |
|---|
| | 249 | // Now decorrelate the equi-variance innovations in preparation for the |
|---|
| | 250 | // log-likelihood computation for this multivariate log-volatility |
|---|
| | 251 | // proposal |
|---|
| | 252 | |
|---|
| 312 | | sum += RVI2(k2, km) * Y2(k2, 6); |
|---|
| 313 | | // ...is the decorrelated multivariate normal that would generate the |
|---|
| 314 | | // next sample's log-volatility given the one proposed for the current |
|---|
| 315 | | // sample. |
|---|
| 316 | | Y2(k2, 7) = sum; |
|---|
| 317 | | } |
|---|
| 318 | | |
|---|
| 319 | | // Then we need to compute the joint log-density of the decorrelated |
|---|
| 320 | | // multivariate normal values |
|---|
| 321 | | sum = - Nz[0] * 0.918938533205; // N*log(sqrt(2*pi)) |
|---|
| | 259 | sum += RI2(k2, km) * Y2(k2, 4); |
|---|
| | 260 | // ...is the decorrelated multivariate innovation for this sample, |
|---|
| | 261 | // given its proposed log-volatility |
|---|
| | 262 | Y2(k2, 5) = sum; |
|---|
| | 263 | } |
|---|
| | 264 | |
|---|
| | 265 | // Now compute log Pr(x[t]|hp[t],h[t+1:]), the log-likelihood of the |
|---|
| | 266 | // decorrelated multivariate innovation, conditional on the proposed |
|---|
| | 267 | // log-volatility hp[t] and the log-volatilities of all subsequent |
|---|
| | 268 | // samples. The proposed log-volatility hp[t] is itself generated |
|---|
| | 269 | // conditional on the log-volatilities of earlier samples, so no further |
|---|
| | 270 | // account of such samples needs be made. |
|---|
| | 271 | |
|---|
| 323 | | for(k2=0; k2<Nz[0]; k2++) |
|---|
| 324 | | sum -= 0.5 * pow(Y2(k2, 7), 2); |
|---|
| 325 | | |
|---|
| 326 | | // Now do the scaling |
|---|
| 327 | | // TODO |
|---|
| 328 | | |
|---|
| 329 | | } |
|---|
| 330 | | |
|---|
| 331 | | |
|---|
| 332 | | |
|---|
| 333 | | |
|---|
| 334 | | |
|---|
| 335 | | |
|---|
| 336 | | |
|---|
| 337 | | |
|---|
| 338 | | |
|---|
| 339 | | // Do a Metropolis-Hastings accept/reject based on the log-likelihood of |
|---|
| 340 | | // the proposal, now in the "sum" variable, versus the log-likelihood of |
|---|
| 341 | | // the previously accepted normal value for this sample. |
|---|
| | 273 | for(k2=0; k2<Nz[0]; k2++) { |
|---|
| | 274 | // Offset the decorrelated innovation by the mean, which is the |
|---|
| | 275 | // correlation between the normal variates for innovation and proposed |
|---|
| | 276 | // volatility shocks, times the value of the normal variate value for |
|---|
| | 277 | // the proposed volatility shock for this sample (or the current normal |
|---|
| | 278 | // variate value, for subsequent samples) |
|---|
| | 279 | if(k3 == k1) { |
|---|
| | 280 | norm = V2(k2, 2); |
|---|
| | 281 | } else { |
|---|
| | 282 | norm = Z2(k2, k3); |
|---|
| | 283 | } |
|---|
| | 284 | Y2(k2, 5) = Y2(k2, 5) - G1(k2) * H2(k2, kh); |
|---|
| | 285 | // The log-likelihood for this sample is simply the argument of the |
|---|
| | 286 | // exponential of the normal pdf, with the reciprocal scaling of the |
|---|
| | 287 | // exponential being done in log space by subtraction. |
|---|
| | 288 | L -= pow(Y2(k2, 5) / Y2(k2, 0), 2) + Y2(k2, 1); |
|---|
| | 289 | // Since the innovation has been inverse-scaled by the square root of |
|---|
| | 290 | // the proposed volatility, the probability density needs to be scaled |
|---|
| | 291 | // as well. In log space, this is a subtraction of half the proposed |
|---|
| | 292 | // log-volatility. |
|---|
| | 293 | L -= 0.5 * H2(k2, kh); |
|---|
| | 294 | } |
|---|
| | 295 | } |
|---|
| | 296 | |
|---|
| | 297 | // TODO |
|---|
| | 298 | // Need to record L for this sample to compare in next MH iteration |
|---|
| | 299 | |
|---|
| | 300 | |
|---|
| | 301 | |
|---|
| | 302 | // Do a Metropolis-Hastings accept/reject based on the log-likelihood for |
|---|
| | 303 | // the proposal versus the log-likelihood for the previously accepted |
|---|
| | 304 | // normal value for this sample. |
|---|
