Changeset 199

Show
Ignore:
Timestamp:
05/27/08 15:37:01 (8 months ago)
Author:
edsuom
Message:

Integration over log-volatility simulations using Metropolis-within-Gibbs is looking most promising

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • projects/AsynCluster/branches/metropolis-within-gibbs/svpmc/model.c

    r198 r199  
    273273// sc   Constants for multivariate normal PDF 
    274274 
    275 //--- Return values (tuple) --- 
    276 // 0    Lx: Log-likelihood of innovations given simulated log-volatilities 
    277 // 1    Lh: Log-likelihood of simulated log-volatilities 
     275//--- Return value (double float) --- 
     276// The linear likelihood P(x|w) of the innovations given the parameters, 
     277// integrated over the ni log-volatility simulations 
    278278 
    279279 
    …… 
    290290int k0, k1, k2, k3; 
    291291double val; 
    292 double Lx = 0; 
    293 double Lh = 0; 
    294292 
    295293// Multivariate innovation for one current time-series sample 
    …… 
    301299double he[2][Nt]; 
    302300// Log probability density of the value of z for current and proposal 
    303 // evaluations 
    304 double Lz[2]; 
     301// evaluations, working values 
     302double Lzk[2]; 
    305303// Log-likelihood of the innovation for the first time-series sample in current 
    306 // and proposal evalutions 
    307 double Lxk[2]; 
     304// and proposal evalutions, working values and accumulated for each iteration 
     305double Lxk[2], Lx[ni]; 
    308306 
    309307// A set of evaluation sample structs 
    …… 
    311309// A struct for the model parameters 
    312310struct parameters P; 
    313  
    314 // The scalar results go here 
    315 py::tuple result(2); 
    316311 
    317312// Generate parameter struct 
    …… 
    328323} 
    329324 
    330  
    331325// For each iteration... 
    332326for(ki=0; ki<ni; ki++) { 
     327 
     328  // Initialize this iteration's accumulator 
     329  Lx[ki] = 0; 
    333330 
    334331  // The "previous" log-volatility of the first time-series sample is set to 
    …… 
    354351    // evaluation struct for the next pair of evaluations. 
    355352    for(k2=0; k2<2; k2++) { 
    356       Lz[k2] = 0; 
    357353      se[k2].Lx = 0; 
     354      Lzk[k2] = 0; 
     355      // Lxk[.] will be set to a value in se[.].Lx that has already been 
     356      // multivariate-accumulated, so it doesn't need clearing here 
    358357    } 
    359358     
    …… 
    379378            // Store the log-density of z for the first time-series sample to 
    380379            // compare current vs. proposal 
    381             Lz[k2] += normLogDensity(val); 
     380            Lzk[k2] += normLogDensity(val); 
    382381          } 
    383382          PA1(se[k2].z, k3) = val; 
    …… 
    411410    // Metropolis-hastings selection of current or proposal using Lz+Lx for 
    412411    // each evaluation 
    413     val = Lz[1] + se[1].Lx  -  Lz[0] - se[0].Lx; 
    414     if(val > 0) { 
    415       // Proposal accepted with probability p=1 
     412    val = Lzk[1] + se[1].Lx  -  Lzk[0] - se[0].Lx; 
     413    if(val > 0 || log(uniform(0, 1)) < val) { 
     414      // Proposal accepted... 
    416415      ke = 1; 
    417     } else if(log(uniform(0, 1)) < val) { 
    418       // Proposal accepted with probability Pr(zp) / Pr(z) 
    419       ke = 1; 
    420       Lh += val; 
     416      // ...overwrite z for this sample with the proposal on z 
     417      for(k1=0; k1<Nt; k1++) 
     418        Z2(k1, k0) = zp[k1]; 
    421419    } else { 
    422420      // Proposal rejected 
    423421      ke = 0; 
    424     } 
    425  
    426     // If the proposal was accepted, overwrite z for this sample with the 
    427     // proposal on z 
    428     if(ke == 1) { 
    429       for(k1=0; k1<Nt; k1++) 
    430         Z2(k1, k0) = zp[k1]; 
    431422    } 
    432423 
    …… 
    440431    } 
    441432 
    442     // Accumulate the log-likelihood of the innovation given the selected 
    443     // log-volatility evaluation 
    444     Lx += Lxk[ke]; 
     433    // The joint probability of the innovation and the variates underlying the 
     434    // log-volatility simulation, conditional on the parameters, is equal to 
     435    // the likelihood of the innovation, conditional on the simulated 
     436    // variates and the parameters, times the probability of the 
     437    // simulated variates, conditional on the parameters: 
     438    // 
     439    // P(x,z|w) = P(x|z,w) * P(z|w) 
     440    // 
     441    // P(z|w) is conditional on w because of the correlation between sequential 
     442    // z values due to the VAR(1) in volatility shocks. Fortunately, the way 
     443    // that we've simulated z allows us to just evaluate P(z) using the joint 
     444    // normal distribution. 
     445 
     446    // Accumulate Log(P(xk,zk|w) = Lxk + Lzk) to eventually get Log(P(x,z|w) 
     447    Lx[ki] += Lxk[ke] + Lzk[ke]; 
    445448   
    446449    // Done with log-volatility draw for this time-series sample 
    …… 
    464467} 
    465468 
    466 // The pair of results is the integrated log-likelihood of the innovations 
    467 // given the log-volatilities simulated over the iterations... 
    468 result[0] = Lx / ni; 
    469 // ...and the integrated log-density of the log-volatilities simulated over the 
    470 // iterations, conditional on the starting z 
    471 result[1] = Lh; 
    472 return_val = result; 
     469// Compute the result, the integrated P(x|w) over the simulations on z 
     470 
     471// Compute the maximum log density // and save it to subtract from everybody 
     472// that the maximum log density is normalized to zero 
     473val = -HUGE_VAL; 
     474for(ki=0; ki<ni; ki++) { 
     475  if(Lx[ki] > val) 
     476    val = Lx[ki]; 
     477} 
     478// Accumulate the linearized normalized densities, borrowing Lzk[0] to do so 
     479Lzk[0] = 0; 
     480for(ki=0; ki<ni; ki++) 
     481  Lzk[0] += exp(Lx[ki] - val); 
     482// ...and return the denormalized, log mean of the linear densities 
     483return_val = log(Lzk[0]) + val - log(ni); 

  • projects/AsynCluster/branches/metropolis-within-gibbs/svpmc/model.py

    r198 r199  
    129129                # Unpack string result from remote call 
    130130                result = list(pack.Unpacker(result)) 
    131             for k, name in enumerate(('Lx', 'Lh', 'z', 'h')): 
     131            for k, name in enumerate(('Lx', 'z', 'h')): 
    132132                setattr(paramContainer, name, result[k]) 
    133133            return paramContainer 
    …… 
    255255        returns a tuple containing the following: 
    256256         
    257             1. The log-likelihood of my observations given the simulated 
    258                log-volatility values, integrated over the simulation rounds. 
    259  
    260             2. The log-probability of the simulated log-volatilities, 
    261                integrated over the simulation rounds. 
    262  
    263             3. The IID normal variates underlying the simulated 
     257            1. The log-likelihood of my observations given the parameters, from 
     258               an integration of the join probability of the observations and 
     259               simulated log-volatilities over the simulation rounds. 
     260                
     261            2. The IID normal variates underlying the simulated 
    264262               log-volatilities after the last simulation round. 
    265263 
    266             4. The multivariate log-volatility value for the last time-series 
     264            3. The multivariate log-volatility value for the last time-series 
    267265               sample after the last simulation round. (Useful for 
    268266               extrapolating from the fitted model.) 
    …… 
    300298 
    301299        # Run the hybrid Gibbs rounds 
    302         Lx, Lh = self.inline( 
     300        Lx = self.inline( 
    303301            'z', 'h', 'x', 'w', 'ni', 
    304302            'd', 'e', 'g', 'rz', 'ri', 'sc', 
    305303            w=self.wiggle, ni=self.Ni) 
    306         return Lx, Lh, z, h 
     304        return Lx, z, h 

  • projects/AsynCluster/branches/metropolis-within-gibbs/svpmc/pmc.py

    r198 r199  
    198198        def weight(paramContainer): 
    199199            if s.isfinite(paramContainer.Lx): 
    200                 L = paramContainer.Lx + paramContainer.Lp \ 
    201                     - paramContainer.Lj - paramContainer.Lh 
     200                L = paramContainer.Lx + paramContainer.Lp - paramContainer.Lj 
    202201                info = " ".join([ 
    203202                    "%+12.2f" % (getattr(paramContainer, x),) 
    204                     for x in ('Lx', 'Lp', 'Lj', 'Lh')]) 
     203                    for x in ('Lx', 'Lp', 'Lj')]) 
    205204            else: 
    206205                info = "" 

  • projects/AsynCluster/branches/metropolis-within-gibbs/svpmc/test/test_model.py

    r198 r199  
    609609        inlineVars.update({'x':2*s.ones((1,Ns)), 'z':s.zeros((1, Ns))}) 
    610610        modelObj = model.Model() 
    611         print "\n  k      Lx       Lh     Simulated Volatility Shocks" 
     611        print "\n  k      Lx       Simulated Volatility Shocks" 
    612612        print "-"*100 
    613613        for k in xrange(Nr): 
    614             Lx, Lh = modelObj.inlineDirect( 
     614            Lx = modelObj.inlineDirect( 
    615615                modelObj.code['likelihood'], **inlineVars) 
    616616            zLine = "   ".join(["%+4.2f" % xx for xx in inlineVars['z'][0,:]]) 
    617             print "%3d: %+6.1f  %+7.1f%s%s" % (k, Lx, Lh, " "*5, zLine) 
     617            print "%3d: %+6.1f  %s%s" % (k, Lx, " "*5, zLine) 
    618618     
    619619    def _simRounds(self, Nr, Ns, inlineVars, name, values): 
    …… 
    636636        """ 
    637637        With univariate, Ni=20 is significantly better than Ni=10, but not 
    638         significantly worse than Ni=40. 
     638        significantly worse than Ni=40. Tested with wiggles of 0.5 and 1.0, and 
     639        eValue=0.95. 
    639640        """ 
    640641        Nr = 1000 
    641         wiggle = 1.0 
     642        wiggle = 0.5 
    642643        eValue = 0.95 
    643644        Ns = 1000 
    …… 
    646647     
    647648    def test_highCorrelation_c(self): 
    648         wiggle = 1.0 
     649        wiggle = 0.5 
    649650        eValue = 0.97 
    650651        Ns, Ni = 200, 20 
    …… 
    654655     
    655656    def test_highCorrelation_model(self): 
    656         wiggle = 1.0 
     657        wiggle = 0.5 
    657658        eValue = 0.97 
    658659        Ns, Ni = 1000, 20 
    …… 
    661662        pc.z = s.zeros((1, Ns)) 
    662663        modelObj.y = self._x_to_y(x) 
    663         Lx, Lh, z, hn = modelObj.likelihood(pc) 
     664        Lx, z, hn = modelObj.likelihood(pc) 
    664665        self._plot(x, (iv[1,:], z[0,:]), (h, self._z_to_h(z[0,:], eValue))) 

  • projects/AsynCluster/trunk/doc/svpmc/example/svpmc.conf

    r196 r199  
    2727Variable    Value           Description 
    2828------------------------------------------------------------------------------- 
    29 D           6               Number of jump deviations in the D-kernel sim 
    30 Df          0.2             Fractional value of each dev from prev (or 1.0) 
     29D           5               Number of jump deviations in the D-kernel sim 
     30Ds          0.1             Starting jump deviation value 
     31Df          0.2             Fractional value of each dev from prev (or start) 
    3132wiggle      0.5             Range +/- of uniform random walks on z for h sim 
    32 N1          30              Proposals on z for each importance-sampling of h 
     33N1          10              Proposals on z for each importance-sampling of h 
    3334N2          20              Hybrid-Gibbs iterations for h sim 
    3435 

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

    r196 r199  
    8686// ---------------------------------------------------------------------------- 
    8787 
     88// Run evaluations until odds of error in selection between proposals null 
     89// proposal on z will be less than 1/100. 
     90#define MIN_DIFF 0.01 
     91 
    8892 
    8993// Model parameters 
    …… 
    107111struct sample 
    108112{ 
     113  // Log-likelihood of the time-series sample in the current evaluation, given 
     114  // h at this point 
     115  double Lxk; 
    109116  // Log-likelihood of the innovations up to the current one, given the history 
    110117  // of h to this point 
    …… 
    223230  // multivariate innovation for this sample, conditional on the sample's 
    224231  // log-volatility. 
    225    
    226   // Now the log-likelihoood. For each time series... 
     232 
     233  S->Lxk = 0; 
     234   
     235  // For each time series... 
    227236  for(k=0; k<N; k++) { 
    228237    // Offset the decorrelated innovation by the mean, which is the correlation 
    …… 
    237246    val = SPA1(S, x0, k) / SPP2(P, sc, k, 0); 
    238247    // Fast squaring instead of pow() 
    239     S->Lx -= 0.5 * val*val + SPP2(P, sc, k, 1); 
     248    S->Lxk -= 0.5 * val*val + SPP2(P, sc, k, 1); 
    240249    // Since the innovation has been inverse-scaled by the square root of the 
    241250    // proposed volatility, the probability density needs to be scaled as 
    242251    // well. In log space, this is a subtraction of half the log-volatility. 
    243     S->Lx -= 0.5 * SPA1(S, h, k); 
    244   } 
    245 } 
    246  
    247  
    248 // Importance-samples one of the proposals based on a weight computed from the 
    249 // log-density of z and the log-likelihood of the innovation given the 
    250 // proposal's log-volatility. 
     252    S->Lxk -= 0.5 * SPA1(S, h, k); 
     253  } 
     254  S->Lx += S->Lxk; 
     255} 
     256 
     257 
     258// Returns with zero only if the non-null proposals have converged close enough 
     259// to the baseline value (the null proposal) to finish likelihood computation 
     260int notCloseEnough(struct sample sp[], int n) 
     261{ 
     262  int k; 
     263  double L; 
     264  for(k=1; k<n; k++) { 
     265    L = sp[k].Lxk - sp[0].Lxk; 
     266    if(L > MIN_DIFF || L < -MIN_DIFF) 
     267      return 1; 
     268  } 
     269  return 0; 
     270} 
     271 
     272 
     273// Importance-samples one of the non-null proposals based on a weight computed 
     274// from the log-density of z and the log-likelihood of the innovation given the 
     275// proposal's log-volatility, conditional on the null proposal. 
    251276// 
    252277// Returns a struct with an integer index to the sampled proposal's struct and 
    …… 
    262287  struct selection result; 
    263288 
    264   // Normalize log-densities to prevent overflow 
    265   for(k=0; k<n; k++) { 
    266     val = sp[k].Lx + Lz[k]; 
     289  // Log-densities densities are relative to the null-proposal, and with a 
     290  // maximum for normalization in selection 
     291  for(k=1; k<n; k++) { 
     292    val = sp[k].Lx + Lz[k] - sp[0].Lx - Lz[0]; 
    267293    if(val > L_max) 
    268294      L_max = val; 
    …… 
    272298  // Compute linear probability density for each proposal, relative to the 
    273299  // maximum one 
    274   for(k=0; k<n; k++) { 
     300  for(k=1; k<n; k++) { 
    275301    val = exp(Lp[k] - L_max); 
    276302    Dp[k] = val; 
    277303    Dp_sum += val; 
    278304  } 
    279  
     305   
    280306  // Accept-reject sampling with uniform random proposal selection 
    281307  do { 
    282308    // Randomly select one of the proposals... 
    283     result.kp = (int) uniform(0, n); 
     309    result.kp = (int) uniform(1, n); 
    284310    // ...until the selected proposal is accepted, with probability 
    285311    // proportional to its relative weight 
    …… 
    305331// ni   Number of hybrid-Gibbs iterations (int) 
    306332// np   Number of proposals to try per log-volatility simulation (int) 
    307 // ne   Max number of time-series samples to evaluate per proposal (int) 
    308333 
    309334//--- Model parameters, direct --- 
    …… 
    385410    // the log-likelihood of the innovation for that sample becomes substantially 
    386411    // the same for all proposals. 
    387    
     412 
    388413    k1 = k0; 
    389414   
    …… 
    405430        for(k3=0; k3<Nt; k3++) { 
    406431          val = Z2(k3, k1); 
    407           if(k1==k0) { 
    408             // First time-series sample, use a proposal on z instead of z itself 
    409             val += uniform(-w, +w); 
     432          if(k1 == k0) { 
     433            // First time-series sample, use a proposal on z instead of z 
     434            // itself... 
     435            if(k2 != 0) 
     436              // ...but only actually change z if this isn't the null proposal 
     437              val += uniform(-w, +w); 
    410438            zp[k2] = val; 
    411439            Lz[k2] = normLogDensity(val); 
    …… 
    432460      k1++; 
    433461     
    434     } while (k1<Ns && k1<k0+ne); 
    435      
     462    } while (k1<Ns && notCloseEnough(sp, np)); 
     463 
    436464    // Importance sample proposal using pProb=exp(Lz+Lx) for each... 
    437465    spSelection = importanceSample(sp, Lz, np); 
    …… 
    480508// simulated log-volatility and Lh, the total log-density of all of the 
    481509// simulated log-volatilities 
    482 result[0] = Lx; 
    483 result[1] = Lh; 
     510result[0] = Lx / ni; 
     511result[1] = Lh / ni; 
    484512return_val = result; 

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

    r196 r199  
    193193      exploration of multivariate space needs more iterations. 
    194194 
    195     @ivar Ne: The number of successive time-series samples to include in the 
    196       evaluation of each log-volatility proposal. For some weird reason, Ne=3 
    197       is optimal, at least for univariate. 
    198      
    199     """ 
    200     keyAttrs = {'y':None, 'wiggle':1.0, 'N1': 30, 'N2':10, 'Ne':3} 
     195    """ 
     196    keyAttrs = {'y':None, 'wiggle':1.0, 'N1': 30, 'N2':10} 
    201197 
    202198    #--- Properties ----------------------------------------------------------- 
    …… 
    309305        sc[:,1] = s.log(s.sqrt(2*s.pi) * sc[:,0]) 
    310306 
    311         print self.wiggle, self.N1, self.N2, self.Ne 
    312307        # Run the hybrid Gibbs rounds, returning the likelihood from the last 
    313308        # one along with the state of the IID variate vector at that point. 
    314309        Lx, Lh = self.inline( 
    315310            'z', 'h', 'x', 
    316             'w', 'ni', 'np', 'ne', 
     311            'w', 'ni', 'np', 
    317312            'd', 'e', 'g', 'rz', 'ri', 'sc', 
    318             w=self.wiggle, np=self.N1, ni=self.N2, ne=self.Ne) 
    319         return Lx/self.N2, Lh/(self.N2*self.Ne), z, h 
     313            w=self.wiggle, np=self.N1+1, ni=self.N2) 
     314        return Lx, Lh, z, h 

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

    r196 r199  
    102102                value = int(value) 
    103103            setattr(self, name, value) 
    104         V = [1.0] 
     104        V = [self.Ds] 
    105105        for k in xrange(self.D-1): 
    106106            V.append(self.Df * V[-1]) 

  • projects/AsynCluster/trunk/svpmc/test/test_model.py

    r192 r199  
    602602        return iv, h, x 
    603603 
    604     def _inlineVars(self, wiggle, eValue, N1, N2=1, Ne=3): 
     604    def _inlineVars(self, wiggle, eValue, N1, N2=1): 
    605605        return { 
    606606            'h':s.empty(1), 
    607607            'sc':s.array([1, s.log(s.sqrt(2*s.pi))]), 
    608608            'w':wiggle, 
    609             'np':N1, 
     609            'np':N1+1, 
    610610            'ni':N2, 
    611             'ne':Ne, 
    612611            'd':s.array([0.0]), 
    613612            'e':s.array([[float(eValue)]]), 
    …… 
    669668        Ns, N1, N2 = 10, 5, 20 
    670669        inlineVars = self._inlineVars(wiggle, eValue, N1, N2) 
    671         inlineVars.update({'x':2*s.ones(Ns), 'z':s.zeros((1, Ns))}) 
     670        inlineVars.update({'x':2*s.ones((1, Ns)), 'z':s.zeros((1, Ns))}) 
    672671        modelObj = model.Model() 
    673672        print "\n  k      Lx       Lh     Simulated Volatility Shocks" 
    …… 
    683682        bins = s.linspace(0.65, 0.95, 20) 
    684683        for k, value in enumerate(values): 
     684            print k, value 
    685685            inlineVars[name] = value 
    686686            correlations = [] 
    …… 
    694694            sp.set_title("%s=%d" % (name, value)) 
    695695 
    696     def test_optimize_Ne(self): 
    697         # Empirically, it looks like Ne=3 is best, just slightly better than 
    698         # Ne=2 and Ne=4. Beyond Ne=4 things decline significantly. 
    699         Nr = 1000 
    700         wiggle = 3.0 
    701         eValue = 0.95 
    702         Ns, N1, N2 = 1000, 10, 2 
    703         inlineVars = self._inlineVars(wiggle, eValue, N1, N2) 
    704         self._simRounds(Nr, Ns, inlineVars, 'ne', (2,3,4,5,6)) 
    705  
    706696    def test_optimize_N1(self): 
    707697        # Empirically, it looks like np=10 is best. Going to np=20 barely gives 
    708698        # any improvement, and obviously doubles the amount of computation 
    709699        # time. Dropping to np=5 is significantly worse. 
    710         Nr = 1000 
    711         wiggle = 3.0 
     700        Nr = 500 
     701        wiggle = 0.5 
    712702        eValue = 0.95 
    713         Ns, N2 = 1000, 2 
     703        Ns, N2 = 500, 10 
    714704        inlineVars = self._inlineVars(wiggle, eValue, 1, N2) 
    715705        self._simRounds(Nr, Ns, inlineVars, 'np', (2,5,10,20)) 
    716706 
    717707    def test_optimize_N2(self): 
    718         # Empirically, it looks like ni=2 is best. It's noticeably better than 
    719         # ni=1 (though not dramatically so), but going to anything larger 
    720         # doesn't show any improvement. 
    721         Nr = 1000 
    722         wiggle = 2.0 
     708        # Empirically, it looks like ni=20 is best. It's noticeably better than 
     709        # ni=10 (though not dramatically so), but going to ni=40 doesn't show 
     710        # any improvement. 
     711        Nr = 500 
     712        wiggle = 0.5 
    723713        eValue = 0.95 
    724         Ns, N1 = 1000, 10 
     714        Ns, N1 = 500, 10 
    725715        inlineVars = self._inlineVars(wiggle, eValue, N1, 1) 
    726         self._simRounds(Nr, Ns, inlineVars, 'ni', (1,2,3,4)) 
     716        self._simRounds(Nr, Ns, inlineVars, 'ni', (5,10,20,40)) 
    727717     
    728718    def test_highCorrelation_c(self): 
    729         Ne = 3 
    730         wiggle = 2.0 
     719        wiggle = 0.5 
    731720        eValue = 0.97 
    732         Ns, N1, N2 = 200, 10, 2 
    733         inlineVars = self._inlineVars(wiggle, eValue, N1, N2, Ne) 
     721        Ns, N1, N2 = 200, 10, 20 
     722        inlineVars = self._inlineVars(wiggle, eValue, N1, N2) 
    734723        for hReal, hSim in self._likelihood(Ns, 1, inlineVars): 
    735724            self._plot((hReal, hSim)) 
    736725     
    737726    def test_highCorrelation_model(self): 
    738         wiggle = 2.0 
    739         eValue = 0.95 
    740         Ns, N1, N2 = 1000, 10, 2 
     727        wiggle = 0.5 
     728        eValue = 0.97 
     729        Ns, N1, N2 = 1000, 10, 20 
    741730        iv, h, x = self._simData(Ns, eValue) 
    742731        pc, modelObj = self._modelAndParamFactory(wiggle, N1, N2, e=[[eValue]])