Changeset 152

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

Working on hybrid Gibbs log-volatility simulation and likelihood computation

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

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

    r147 r152  
    8383 
    8484 
    85 // Model.draw_h 
     85// Model.likelihood 
    8686//  
    8787// Supplied variables 
    8888// ---------------------------------------------------------------------------- 
    89 // h    Volatility values (values to be overwritten), [pn] array 
    90 // v    Volatility shock values, [pn] array 
     89// L    Log-likelihood of innovations, given params & simulated volatilities 
     90 
     91//--- Supplied arrays --- 
     92// z    Independent normal variates, [pn] array 
     93// x    Innovation values, [pn] array 
     94// r    Concurrent cross-correlations between volatility shocks, [pp] 
     95// s    Inverse, concurrent cross-correlations between innovation shocks, [pp] 
     96 
     97//--- Externally declared scratchpad --- 
     98// v    Volatility shock values, [p] vector for one multivariate sample 
     99// h    Current and previous log-volatilities, [p2] array 
     100// y    Volatility-normalized and decorrelated innovation values, [p2] array 
     101 
     102//--- Model parameters --- 
    91103// d    Volatility offset, [p] vector 
    92 // e    Lag-1 cross-correlation [pp] array 
    93  
    94 int i, j, k; 
    95 double sum; 
    96   
    97 // The first volatility value for each time series is just the first volatility 
    98 // shock plus the offset 
    99 for(i=0; i<Nv[0]; i++) 
    100   H2(i,0) = V2(i,0) + D1(i); 
    101 // For each innovation after the first... 
    102 for(j=1; j<Nv[1]; j++) { 
    103   // For each time series... 
    104   for(i=0; i<Nv[0]; i++) { 
    105     // The offset plus the current shock... 
    106     sum = D1(i) + V2(i,j); 
    107     // ...plus the dot product for the VAR(1) term... 
    108     for(k=0; k<Nv[0]; k++) 
    109       sum += E2(i,k) * H2(k,j-1); 
    110     // ...is the modeled output 
    111     H2(i,j) = sum; 
     104// e    Lag-1 cross-correlations for VAR of volatilities, [pp] array 
     105// g    Volatility/innovation shock correlations, [p] vector 
     106 
     107//--- Scalar parameters --- 
     108// n0   Number of volatility draws per likelihood evaluation.   
     109// n3   Number of random-walk normal draws per volatility draw. 
     110// wiggle = Uniform proposal +/- extent 
     111 
     112 
     113// THIS IS WHERE THE VAST MAJORITY OF THE CPU TIME IS EXPENDED! 
     114 
     115 
     116int km, kh; 
     117int k0, k1, k2, k3; 
     118double sum, normp, dL, L_normp, Lp; 
     119 
     120 
     121// For each overall iteration... 
     122for(k0=0; k0<n0; k0++) { 
     123 
     124  // The log-likelihood of the innovations given the log-volatility proposal to 
     125  // be constructed 
     126  Lp = 0; 
     127 
     128  // Initialize log-volatility scratchpad: the "previous" volatility value 
     129  // for the first multivariate sample is just the offset 
     130  for(k2=0; k2<Nz[0]; k2++) 
     131      H2(k2, 1) = D1(i); 
     132 
     133  // for each multivariate sample... 
     134  for(k1=0; k1<Nz[1]; k1++) { 
     135 
     136    // Update this column of independent random-walk normal variates to produce 
     137    // a multivariate i.i.d. normal proposal. It is this on which we build a 
     138    // multivariate log-volatility proposal for this multivariate sample. 
     139 
     140    // For each time series... 
     141    for(k2=0; k2<Nz[0]; k2++) { 
     142      // Get the log-likelihood of the current value 
     143      L_norm = -0.5 * pow(Z2(k2, k1), 2); 
     144      // For each oversampling iteration... 
     145      for(k3=0; k3<n3; k3++) { 
     146        // Generate a uniform, symmetric, random walk proposal 
     147        normp = Z2(k2, k1) + wiggle * (drand48() - 0.5); 
     148        // Do the ratio in log space 
     149        L_normp = -0.5 * pow(normp, 2); 
     150        dL = L_normp - L_norm; 
     151        // The uniform random variate thus needs to be in logspace as well, but 
     152        // it's not always even needed. Thus the log operation is done just 
     153        // once, and only some of the time. 
     154        if(dL > 0  || log(drand48()) < dL) { 
     155          // Update current value and its log probability when proposal 
     156          // accepted 
     157          Z2(k2, k1) = normp; 
     158          L_norm = L_normp; 
     159        } 
     160      } 
     161    } 
     162 
     163    // Now correlate the i.i.d. proposal to form a proposed column of 
     164    // volatility shocks 
     165 
     166    // For each time series... 
     167    for(k2=0; k2<Nz[0]; k2++) { 
     168      // The matrix product of the cross-correlations and the independent 
     169      // variates... 
     170      sum = 0; 
     171      for(km=0; km<Nz[0]; km++) 
     172        sum += R2(k2, km) * Z2(k2, k1); 
     173      // ...is the correlated volatility shock for this time series 
     174      V1(k2) = sum; 
     175    } 
     176 
     177    // Now do the VAR(1) to get a multivariate log-volatility proposal 
     178 
     179    // We alternate between the two columns of the log-volatility scratchpad 
     180    // for current & previous multivariate log-volatility values 
     181    kh = k1 % 2; 
     182 
     183    // For each time series... 
     184    for(k2=0; k2<Nv[0]; k2++) { 
     185      // The offset plus the current shock... 
     186      sum = D1(i) + V1(k2); 
     187      // ...plus the dot product for the VAR(1) term... 
     188      for(km=0; km<Nv[0]; km++) 
     189        sum += E2(k2, km) * H2(k2, 1-kh); 
     190      // ...is the modeled output 
     191      H2(k2, kh) = sum; 
     192    } 
     193 
     194    // Now inverse-scale the innovations for this sample by the 
     195    // log-volatilities in preparation for decorrelating them 
     196    for(k2=0; k2<Nv[0]; k2++) 
     197      Y2(k2, 0) = exp(-0.5 * H2(k2, kh)) * Z2(k2, k1); 
     198     
     199    // Now decorrelate the equi-variance innovations in preparation for 
     200    // (finally) the log-likelihood computation for this multivariate sample 
     201     
     202    // For each time series... 
     203    for(k2=0; k2<Nz[0]; k2++) { 
     204      // The product of the inverted cross-correlation matrix and the 
     205      // correlated variates... 
     206      sum = 0; 
     207      for(km=0; km<Nz[0]; km++) 
     208        sum += S2(k2, km) * Y2(k2, 0); 
     209      // ...is the deccorrelated multivariate innovation for this sample 
     210      Y2(k2, 1) = sum; 
     211    } 
     212 
     213    // Finally, compute the log-likelihood of this deccorrelated multivariate 
     214    // innovation 
     215 
     216    // TODO 
     217     
    112218  } 
    113219} 
    114  
    115  
    116  
    117 // Model.decorrelate 
    118 // 
    119 // Supplied variables 
    120 // ---------------------------------------------------------------------------- 
    121 // y    2-D array of correlated random variates (to be overwritten) 
    122 // x    2-D array of independent random variates 
    123 // r    2-D array of inverted correlations between values in each column of x 
    124  
    125 int i, j, k; 
    126 double sum; 
    127  
    128 // For each column of independent variates to be decorrelated... 
    129 for(i=0; i<Nx[1]; i++) { 
    130   // For each of the values in the column... 
    131   for(j=0; j<Nx[0]; j++) { 
    132     // The product of the inverted cross-correlation matrix and the correlated 
    133     // variates... 
    134     sum = 0; 
    135     for(k=0; k<Nx[0]; k++) 
    136       sum += R2(j,k) * X2(k,i); 
    137     // ...is the correlated output 
    138     Y2(j, i) = sum; 
    139   } 
    140 } 

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

    r151 r152  
    180180      observation data. 
    181181 
    182     @ivar M: Number of volatility draws per likelihood evaluation. 
     182    @ivar n0: Number of volatility draws per likelihood evaluation. 
     183 
     184    @ivar n3: Number of random-walk normal draws per volatility draw. 
    183185 
    184186    """ 
    185187    _logU_index = -1 
    186188 
    187     keyAttrs = {'tsList':None, 'M':500} 
     189    keyAttrs = {'tsList':None, 'n0':500, 'n3':10} 
    188190    paramNames = params.PARAM_NAMES 
    189191 
    …… 
    278280        return logProbs.sum() 
    279281 
    280     def logUniform(self): 
    281         """ 
    282         Returns a log-uniform variate taken from a large array that is 
    283         generated and refreshed periodically. Generating the uniform variates 
    284         and taking their log in a single large array operation is more 
    285         efficient than doing so one value at a time. 
    286         """ 
    287         logU = getattr(self, '_logU_array', []) 
    288         if self._logU_index >= len(logU)-1: 
    289             self._logU_index = -1 
    290             # The efficient array operation 
    291             logU = self._logU_array = s.log(s.rand(10000)) 
    292         self._logU_index += 1 
    293         return logU[self._logU_index] 
     282    def covarMatrix(self, correlations): 
     283        """ 
     284        Returns a covariance matrix from the vector or sequence of 
     285        I{correlations} supplied in the form used by L{correlate}. 
     286 
     287        Assumes that the independent components have unit variance. 
     288        """ 
     289        i = 0 
     290        p = self.p 
     291        cv = s.ones((p, p)) 
     292        for j in xrange(p - 1): 
     293            for k in xrange(j+1, p): 
     294                cv[j, k] = cv[k, j] = correlations[i] 
     295                i += 1 
     296        return cv 
    294297 
    295298 
    296299    #--- The API -------------------------------------------------------------- 
    297300 
    298     def likelihood(self, paramContainer, wiggle): 
     301    def likelihood(self, paramContainer, sigma): 
    299302        """ 
    300303        Call this method with a parameter object that defines: 
    …… 
    302305            1. a set of values for my parameters, and 
    303306 
    304             2. a latent parameter consisting of an array of volatility shocks 
    305               to be used as a starting point for another simulation draw of 
    306               volatilities, and 
     307            2. a latent parameter consisting of an array of i.i.d. normal 
     308              variates to be used as a starting point for another simulation 
     309              draw of log volatilities, and 
    307310 
    308311            3. a 1-D array of log-likelihoods for each observation time given 
    309                the volatility shocks. 
    310  
    311         A new 2-D array of log-volatilities will be drawn from a relatively 
    312         short Metropolis-Hastings random walk simulation of the volatility 
    313         shocks coupled with a VAR(1) from the shocks to the 
    314         log-volatilities. The Metropolis-Hastings proposals are scaled by the 
    315         supplied fractional I{wiggle}. The simulation is governed by the 
     312               the i.i.d. variates after cross-correlation and VAR(1). 
     313 
     314        A new 2-D array of log-volatilities will be drawn from a hybrid Gibbs 
     315        simulation, with each Gibbs step driven by a short Metropolis-Hastings 
     316        random walk simulation of the i.i.d. variates, followed by 
     317        cross-correlation to form volatility shocks, and VAR(1) to form the 
     318        log-volatilities. The Metropolis-Hastings proposals are scaled based on 
     319        the supplied fractional I{wiggle}. The simulation is governed by the 
    316320        log-likelihood of my observation data given the supplied parameters and 
    317         the simulated volatilities. 
    318  
    319         The value of that log-likelihood at the end of the log-volatility 
    320         simulation is returned in a tuple, followed by the volatility shocks 
    321         underlying the log-volatilities drawn at that point. 
    322  
     321        the simulated log-volatilities. 
     322 
     323        The value of the log-likelihood at the end of the log-volatility 
     324        simulation is returned in a tuple, followed by the i.i.d. normal 
     325        variates underlying the log-volatilities drawn at that point. 
     326         
    323327        The returned likelihood does not consider the prior probability of the 
    324328        parameters. You have to account for that yourself, preferably before 
    …… 
    327331        what, making any call to this method with it a waste of computing time. 
    328332        """ 
    329         if len(paramContainer.latent): 
    330             self.nw(paramContainer.latent) 
     333        # Set my parameters... 
    331334        self.setParamSequence(paramContainer.paramSequence()) 
    332         # Obtain innovations as residuals of the reversed VAR(1) 
     335        # ...including the latent parameter of i.i.d. normal variates 
     336        z = paramContainer.latent 
     337        # Innovations as residuals of the reversed VAR(1) 
    333338        x = self.var.reverse(self.a, self.b) 
    334         # MCMC simulation of log-volatilities 
    335         L = -s.inf 
    336         for j in xrange(self.M): 
    337             # Draw a proposal array of log-volatilities 
    338             h = self.draw_h(wiggle) 
    339             # Compute the total log-likelihood conditional on the proposed 
    340             # log-volatilities 
    341             #Lp = self.likelihoodOfInnovations(s.exp(-0.5*h)*x) 
    342             Lp = self.likelihoodOfInnovations(self.decorrelate(x, h)) 
    343             # Metropolis-Hastings accept/reject of the proposal 
    344             if Lp > L or (Lp-L) > self.logUniform(): 
    345                 if j > 0: 
    346                     print "\n%d -> %d" % (L, Lp) 
    347                 L = Lp 
    348                 vShocks = self.nw() 
    349             else: 
    350                 print ".", 
    351         print "\n" 
    352         return L, vShocks 
     339        # Concurrent cross-correlations between volatility shocks 
     340        r = linalg.cholesky(self.covarMatrix(self.f), lower=True)) 
     341        # Inverse of concurrent cross-correlations between innovation shocks 
     342        s = linalg.inv(linalg.cholesky(self.covarMatrix(self.c), lower=True)) 
     343        # Declare some scratchpad arrays 
     344        v = s.empty(self.p) 
     345        h = s.empty((self.p, 2)) 
     346        y = s.empty((self.p, 2)) 
     347        # Finally, run the hybrid Gibbs simulation of log-volatilities to get a 
     348        # log-likelihood that is not (very) conditional on the latent 
     349        # parameters that the volatilities represent 
     350        L = self.inline( 
     351            'L', 
     352            'v', 'x', 'r', 's', 
     353            'z', 'h', 'y', 
     354            'd', 'e', 'g', 
     355            'n0', 'n3', 'wiggle', 
     356            L=0, wiggle=sigma/s.sqrt(self.n3)) 
     357        return L, v 
     358     
     359