Changeset 155

Show
Ignore:
Timestamp:
04/23/08 02:17:19 (8 months ago)
Author:
edsuom
Message:

Got hybrid Gibbs Model.likelihood prototyped

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • projects/AsynCluster/trunk/svpmc/model.c

    r154 r155  
    8383 
    8484 
    85 // Model.likelihood 
     85// Model.zProposal 
     86// 
     87// Supplied variables 
     88// ---------------------------------------------------------------------------- 
     89// z    Independent normal variates, [pn] array (updated) 
     90// M    Number of uniform-proposal MCMC iterations per IID normal draw 
     91// wiggle = Width of random-walk uniform proposal 
     92 
     93 
     94int k0, k1, k2; 
     95double norm, normp, L, Lp; 
     96 
     97 
     98// for each multivariate sample... 
     99for(k0=0; k0<Nz[1]; k0++) { 
     100 
     101  // Do a random walk from this column of independent random-walk normal 
     102  // variates to produce a multivariate iid normal proposal. It is this on 
     103  // which we build a multivariate log-volatility proposal for this 
     104  // multivariate sample. 
     105 
     106  // For each time series... 
     107  for(k1=0; k1<Nz[0]; k1++) { 
     108    // Initialize the current normal value with the current proposal value and 
     109    // compute its log-likelihood 
     110    norm = Z2(k1, k0); 
     111    L = -0.5 * pow(norm, 2); 
     112    // For each oversampling iteration... 
     113    for(k2=0; k2<M; k2++) { 
     114      // Generate a uniform, symmetric, random walk proposal 
     115      normp = norm + wiggle * (drand48() - 0.5); 
     116      // Do the ratio in log space 
     117      Lp = -0.5 * pow(normp, 2); 
     118      dL = Lp - L; 
     119      // The uniform random variate thus needs to be in logspace as well, but 
     120      // it's not always even needed. Thus the log operation is done just 
     121      // once, and only some of the time. 
     122      if(dL > 0  || log(drand48()) < dL) { 
     123        // Update current value and its log probability when proposal 
     124        // accepted 
     125        norm = normp; 
     126        L = Lp; 
     127      } 
     128    } 
     129    // Update this element of the IID normal array in place 
     130    Z2(k1, k0) = norm; 
     131  } 
     132} 
     133 
     134 
     135 
     136// Model.logVolatilities 
    86137//  
    87138// Supplied variables 
    …… 
    89140 
    90141//--- Supplied arrays --- 
    91 // z    Independent normal variates for log-volatilities, [pn] array (updated) 
    92 // v    Log-volatility shocks, [pn] (overwritten) 
     142// z    Independent normal variates for log-volatilities, [pn] array 
     143// v    Log-volatility shocks, [pn] 
    93144// x    Innovation values, [pn] array 
    94 // rv   Concurrent cross-correlations between volatility shocks, [pp] 
     145// h    Log-volatility values, [p(n+1)] array (last n samples overwritten) 
    95146// ri   Inverse, concurrent cross-correlations between innovation shocks, [pp] 
    96147// y    Miscellaneous multivariate scratchpad, [p8] array (overwritten) 
    97 // h    Log-volatility scratchpad, [p2] array, for VAR(1) 
    98148 
    99149//--- Model parameters --- 
    …… 
    102152// g    Volatility/innovation shock correlations, [p] vector 
    103153 
    104 //--- Scalar parameters --- 
    105 // n0   Number of volatility draws per likelihood evaluation.   
    106 // n3   Number of random-walk normal draws per volatility draw. 
    107 // nv   Number of samples after current one to compute VAR effect 
    108 // sigma = Standard deviation of random-walk normal proposal 
    109  
    110154 
    111155// The multivariate scratchpad y is used as follows, by column: 
    …… 
    113157// 0    2*(1-sigma**2) for bivariate normal probability density function 
    114158// 1    0.5*log(pi*2*sigma**2) offset for bivariate normal pdf 
    115 // 2    Independent normal variate random-walk proposals 
    116 // 3    Correlated volatility shocks 
    117 // 4    Innovation values after inverse-scaling by log-volatility 
    118 // 5    Inverse-scaled innovation values after decorrelation 
    119 // 6    - 
    120 // 7    - 
    121  
    122  
    123 // THIS CODE IS WHERE THE VAST MAJORITY OF THE TOTAL CPU TIME IS EXPENDED! 
    124  
    125  
    126 int km, kh; 
    127 int k0, k1, k2, k3, nd; 
    128 double sum, norm, normp, dL, L0, L1, L; 
    129  
    130 double wiggle = sigma * pow(n3, -0.5); 
    131  
    132  
    133 // Initialize  constants for bivariate normal probability density function 
    134 for(k2=0; k2<Nz[0]; k2++) { 
    135   Y2(k2, 0) = 2*(1 - pow(G1(k2), 2)); 
    136   Y2(k2, 1) = 0.5 * log(3.1415926535897931 * Y2(k2, 0)); 
    137 } 
    138  
    139  
    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... 
     159// 2    Innovation values after inverse-scaling by log-volatility 
     160// 3    Inverse-scaled innovation values after decorrelation 
     161 
     162 
     163int k0, k1, km; 
     164double sum; 
     165double L = 0; 
     166 
     167 
     168// For each sample... 
     169for(k0=0; k0<Nz[1]; k0++) { 
     170 
     171  // Run the VAR(1) process to get the current multivariate log-volatility 
     172  // value 
     173 
     174  // For each time series... 
     175  for(k1=0; k1<Nz[0]; k1++) { 
     176    // The volatility offset plus the shock... 
     177    sum = D1(k1) + V2(k2, k1); 
     178    // ...plus the dot product for the VAR(1) term... 
     179    for(km=0; km<Nz[0]; km++) 
     180      sum += E2(k1, km) * H2(k1, k0); 
     181    // ...is the current log-volatility value 
     182    H2(k1, k0+1) = sum; 
     183  } 
     184     
     185  // Now inverse-scale the innovations for this sample by the square root of 
     186  // the volatilities in preparation for decorrelating the innovations 
     187  for(k1=0; k1<Nz[0]; k1++) 
     188    Y2(k1, 2) = exp(-0.5 * H2(k1, k0+1)) * X2(k1, k0); 
     189     
     190  // Now decorrelate the equi-variance innovations in preparation for the 
     191  // log-likelihood computation for this multivariate log-volatility proposal 
     192     
     193  // For each time series... 
     194  for(k1=0; k1<Nz[0]; k1++) { 
     195    // The product of the inverted cross-correlation matrix and the correlated 
     196    // variates... 
    150197    sum = 0; 
    151198    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  
    162 // For each overall iteration... 
    163 for(k0=0; k0<n0; k0++) { 
    164  
    165   // for each multivariate sample... 
    166   for(k1=0; k1<Nz[1]; k1++) { 
    167  
    168     // Do a random walk from this column of independent random-walk normal 
    169     // variates to produce a multivariate iid normal proposal. It is this on 
    170     // which we build a multivariate log-volatility proposal for this 
    171     // multivariate sample. 
    172  
    173     // For each time series... 
    174     for(k2=0; k2<Nz[0]; k2++) { 
    175       // Initialize the scratchpad with the current proposal value and compute 
    176       // its log-likelihood 
    177       Y2(k2, 2) = Z2(k2, k1); 
    178       L0 = -0.5 * pow(Y2(k2, 2), 2); 
    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; 
    194         } 
    195       } 
    196     } 
    197      
    198     // Now correlate the iid proposal to form a proposed column of volatility 
    199     // shocks 
    200      
    201     // For each time series... 
    202     for(k2=0; k2<Nz[0]; k2++) { 
    203       // The matrix product of the cross-correlations and the proposed 
    204       // independent variates... 
    205       sum = 0; 
    206       for(km=0; km<Nz[0]; km++) 
    207         sum += RV2(k2, km) * Y2(k2, 2); 
    208       // ...is the correlated volatility shock proposal for this time series 
    209       Y2(k2, 3) = sum; 
    210     } 
    211      
    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; 
    220     } 
    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; 
    227       // For each time series... 
    228       for(k2=0; k2<Nz[0]; k2++) { 
    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      
    253       // For each time series... 
    254       for(k2=0; k2<Nz[0]; k2++) { 
    255         // The product of the inverted cross-correlation matrix and the 
    256         // correlated variates... 
    257         sum = 0; 
    258         for(km=0; km<Nz[0]; km++) 
    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. 
     199      sum += RI2(k1, km) * Y2(k1, 2); 
     200    // ...is the decorrelated multivariate innovation for this sample, given 
     201    // its log-volatility 
     202    Y2(k1, 3) = sum; 
     203  } 
     204     
     205  // Now compute log Pr(x[t]|h[t]), the log-likelihood of the decorrelated 
     206  // multivariate innovation for this sample, conditional on the just-computed 
     207  // log-volatility h[t]. 
    271208       
    272       // For each time series... 
    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. 
    305  
    306     // If this is the first, initialization iteration, we always accept the 
    307     // null "proposal," which was only obtained to (re)compute the 
    308     // log-likelihood given the current normal variates. 
    309     if(k0 == 0) { 
    310       P2(0, k1) = L0; 
    311       P2(1, k1) = L1; 
    312     } else { 
    313       // Do the ratio of likelihoods for hp[t] vs current h[t] in log space 
    314       dL = L1 - P2(0, k1); 
    315       // The uniform random variate thus needs to be in logspace as well, but 
    316       // it's not always even needed. Thus the log operation is done just once, 
    317       // and only some of the time. 
    318       if(dL > 0 || log(drand48()) < dL) { 
    319         P2(0, k1) = L0; 
    320         P2(1, k1) = L1; 
    321         // For each time series... 
    322         for(k2=0; k2<Nz[0]; k2++) 
    323           Z2(k2, k1) = Y2(k2, 2); 
    324       } 
    325     } 
    326  
    327     // End of Metropolis-Hastings step for this multivariate sample 
    328   } 
    329  
    330   // End of this iteration 
    331 } 
    332  
    333 // End of MCMC run. The array Z2 has been overwritten with updated iid normal 
    334 // variates underlying an updated set of log-volatilities, and the vector P1 
    335 // has been overwritten with the log-likelihood of each multivariate innovation 
    336 // sample given those log-volatilities. 
    337  
    338  
     209  // For each time series... 
     210  for(k1=0; k1<Nz[0]; k1++) { 
     211    // Offset the decorrelated innovation by the mean, which is the correlation 
     212    // between the normal variates for innovation and proposed volatility 
     213    // shocks, times the value of the normal variate value underlying the 
     214    // volatility shock for this sample 
     215    Y2(k1, 3) = Y2(k1, 3)  -  G1(k1) * Z2(k1, k0); 
     216    // The log-likelihood for this sample is simply the argument of the 
     217    // exponential of the normal pdf, with the reciprocal scaling of the 
     218    // exponential being done in log space by subtraction. 
     219    L -= pow(Y2(k1, 3) / Y2(k1, 0), 2) + Y2(k1, 1); 
     220    // Since the innovation has been inverse-scaled by the square root of the 
     221    // proposed volatility, the probability density needs to be scaled as 
     222    // well. In log space, this is a subtraction of half the log-volatility. 
     223    L -= 0.5 * H2(k1, k0+1); 
     224  } 
     225} 
     226 
     227// The return value is the total log-volatility 
     228return_val = L; 

  • projects/AsynCluster/trunk/svpmc/model.py

    r154 r155  
    184184      observation data. 
    185185 
    186     @ivar n0: Number of volatility draws per likelihood evaluation. 
    187  
    188     @ivar n3: Number of random-walk normal draws per volatility draw. 
     186    @ivar Mv: Number of volatility draws per likelihood evaluation. 
     187 
     188    @ivar Mz: Number of random-walk normal draws per volatility draw. 
    189189 
    190190    """ 
    …… 
    226226    #--- Support -------------------------------------------------------------- 
    227227 
    228     def decaySamples(self, tol): 
    229         """ 
    230         """ 
    231          
    232  
    233  
    234228    def covarMatrix(self, correlations): 
    235229        """ 
    …… 
    248242        return cv 
    249243 
     244    def decaySamples(self, tol): 
     245        """ 
     246        """ 
     247 
     248    def logUniform(self): 
     249        """ 
     250        Returns a log-uniform variate taken from a large array that is 
     251        generated and refreshed periodically. Generating the uniform variates 
     252        and taking their log in a single large array operation is more 
     253        efficient than doing so one value at a time. 
     254        """ 
     255        logU = getattr(self, '_logU_array', []) 
     256        if self._logU_index >= len(logU)-1: 
     257            self._logU_index = -1 
     258            # The efficient array operation 
     259            logU = self._logU_array = s.log(s.rand(10000)) 
     260        self._logU_index += 1 
     261        return logU[self._logU_index] 
     262 
     263    def logVolatilities(self, z, v, x, h0): 
     264        """ 
     265        Call this method with a 2-D array I{z} of multivariate IID normal 
     266        random variates, a 2-D array I{v} of volatility shocks correlated from 
     267        the IID variates, a 2-D vector I{x} containing one multivariate 
     268        innovation for each random variate in I{z}, and a 1-D vector I{h0} that 
     269        defines an initial multivariate log-volatility value. 
     270 
     271        The method will compute a 2-D array I{h} of log-volatilities for each 
     272        volatility shock by running it through the stochastic-volatility VAR(1) 
     273        process, with the value provided in I{h0} as the VAR component for the 
     274        first sample. Then it will compute the total log-likelihood of the 
     275        innovations in I{x} given each respective log-volatility value in I{h}. 
     276         
     277        The method returns a tuple with the total log-likelihood value and the 
     278        computed log-volatility array I{h}. 
     279        """ 
     280        # Mostly empty log-volatility array 
     281        h = s.column_stack([h0, s.empty_like(v)]) 
     282        # Working arrays with intialized values 
     283        if 'lv_work' in self.cache: 
     284            y, ri = self.cache['lv_work'] 
     285        else: 
     286            # Scratchpad with some constants 
     287            y = s.empty((self.p, 4)) 
     288            y[:,0] = 2 * (1 - self.g**2) 
     289            y[:,1] = 0.5 * s.log(s.pi*y[:,0]) 
     290            # Inverse of concurrent cross-correlations between innovation shocks 
     291            ri = linalg.inv( 
     292                linalg.cholesky(self.covarMatrix(self.c), lower=True)) 
     293            self.cache['lv_work'] = y, ri 
     294        # Run the C code where the heavy lifting is done 
     295        L = self.inline( 
     296            'z', 'v', 'x', 'h', 
     297            'ri', 'y', 'd', 'e', 'g') 
     298        return L, h[:,1:] 
     299 
     300    def zProposal(self, z, sigma): 
     301        """ 
     302        """ 
     303        z = z.copy() 
     304        self.inline( 
     305            'z', 'M', 'wiggle', M=self.Mz, wiggle=sigma/s.sqrt(self.Mz)) 
     306        return z 
     307     
     308    def hybridGibbs(self, z, x, sigma): 
     309        """ 
     310        Does one iteration of a hybrid Gibbs sampler for the log-volatilities, 
     311        where each conditional draw is done with a Metropolis-Hastings step. 
     312         
     313        Returns the likelihood of the innovations in I{x} given the simulated 
     314        log-volatilities. Updates the IID normal variates in place. 
     315        """ 
     316        L_total = 0 
     317        if 'hg_work' in self.cache: 
     318            rv, Nv = self.cache['hg_work'] 
     319        else: 
     320            # Concurrent cross-correlations between volatility shocks 
     321            rv = s.matrix( 
     322                linalg.cholesky(self.covarMatrix(self.f), lower=True)) 
     323            Nv = 100 # TODO: Use decaySamples() 
     324            self.cache['hg_work'] = rv, Nv 
     325        # Initialize volatility shocks from cross-correlation of the IID 
     326        # variates... 
     327        v = rv * z 
     328        # ...and initial log-volatilites, from the offset alone... 
     329        h0 = self.d.copy() 
     330        # ...and the log-likelihood of the first sample given the current IID 
     331        # variates... 
     332        L = self.logVolatilities(z[:,:N], v[:,:N], x[:,:N], h0)[0] 
     333        # ...and a set of random-walk proposals for the IID variates, along 
     334        # with their volatility shocks 
     335        zp = self.zProposal(z, sigma) 
     336        vp = rv * zp; 
     337        # Do conditional draws of each sample in turn 
     338        for k in xrange(self.n): 
     339            zpk = s.column_stack([zp[:,k], z[:,k+1:k+N]]) 
     340            vpk = s.column_stack([vp[:,k], v[:,k+1:k+N]]) 
     341            Lp, h = self.logVolatilities(zpk, vpk, x[:,k:k+N], h0) 
     342            # Metropolis-Hastings accept/reject of the proposal 
     343            if Lp > L or (Lp-L) > self.logUniform(): 
     344                z[k] = zp[k] 
     345                v[k] = vp[k] 
     346                h0 = h[0] 
     347            # Add the log-likelihood of the drawn log-volatility sample to the 
     348            # overall tally 
     349            L_total += L 
     350        return L_total 
     351     
    250352 
    251353    #--- The API -------------------------------------------------------------- 
    252  
     354     
    253355    def likelihood(self, paramContainer, sigma): 
    254356        """ 
    …… 
    257359            1. a set of values for my parameters, and 
    258360 
    259             2. a latent parameter consisting of an array of i.i.d. normal 
     361            2. a latent parameter consisting of an array of IID normal 
    260362              variates to be used as a starting point for another simulation 
    261363              draw of log volatilities, and 
    262364 
    263365            3. a 1-D array of log-likelihoods for each observation time given 
    264                the i.i.d. variates after cross-correlation and VAR(1). 
     366               the IID variates after cross-correlation and VAR(1). 
    265367 
    266368        A new 2-D array of log-volatilities will be drawn from a hybrid Gibbs 
    267369        simulation, with each Gibbs step driven by a short Metropolis-Hastings 
    268         random walk simulation of the i.i.d. variates, followed by 
     370        random walk simulation of the IID variates, followed by 
    269371        cross-correlation to form volatility shocks, and VAR(1) to form the 
    270372        log-volatilities. The Metropolis-Hastings proposals are scaled based on 
    …… 
    274376 
    275377        The value of the log-likelihood at the end of the log-volatility 
    276         simulation is returned in a tuple, followed by the i.i.d. normal 
     378        simulation is returned in a tuple, followed by the IID normal 
    277379        variates underlying the log-volatilities drawn at that point. 
    278380         
    …… 
    285387        # Set my parameters... 
    286388        self.setParamSequence(paramContainer.paramSequence()) 
    287         # ...including the latent parameter of i.i.d. normal variates 
    288         z = paramContainer.latent 
     389        # ...including the latent parameter of IID normal variates 
     390        z = paramContainer.latent.copy() 
    289391        # Innovations as residuals of the reversed VAR(1) 
    290392        x = self.var.reverse(self.a, self.b) 
    291         # Concurrent cross-correlations between volatility shocks 
    292         rv = linalg.cholesky(self.covarMatrix(self.f), lower=True) 
    293         # Inverse of concurrent cross-correlations between innovation shocks 
    294         ri = linalg.inv(linalg.cholesky(self.covarMatrix(self.c), lower=True)) 
    295         # Declare working array with zeros 
    296         y = s.zeros((self.p, 8)) 
    297         # Finally, run the hybrid Gibbs simulation of log-volatilities to get a 
    298         # log-likelihood that is not (very) conditional on the latent 
    299         # parameters that the volatilities represent 
    300         v = s.empty(*self.p, self.n) 
    301         h = s.empty((self.p, 2)) 
    302         L = self.inline( 
    303             'z', 'v', 'x', 'rv', 'ri', 'y', 'h', 
    304             'd', 'e', 'g', 
    305             'n0', 'n3', 'sigma') 
     393        for k in xrange(self.Mv): 
     394            L = self.hybridGibbs(z, x, sigma) 
    306395        return L, z 
    307      
    308     def foo(self): 
    309         pass