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

/* svPMC: Stochastic Volatility Inference via Population Monte Carlo
   
   Copyright (C) 2007-2008 by Edwin A. Suominen, http://www.eepatents.com
   
   This program is free software; you can redistribute it and/or modify it under
   the terms of the GNU General Public License as published by the Free Software
   Foundation; either version 2 of the License, or (at your option) any later
   version.
   
   This program is distributed in the hope that it will be useful, but WITHOUT
   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
   FOR A PARTICULAR PURPOSE.  See the file COPYING for more details.
   
   You should have received a copy of the GNU General Public License along with
   this program; if not, write to the Free Software Foundation, Inc., 51
   Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA

*/


// VAR.reverse
// 
// Supplied variables
// ----------------------------------------------------------------------------
// v    Innovation values (initially empty), [pn] array
// y    Observation values, [pn] array
// a    Drift, [p] vector
// b    Lag-1 cross-correlation [pp] array

int i, j, k;
double sum;

// The first innovation for each time series is just the first observation
// minus the drift
for(i=0; i<Ny[0]; i++)
  V2(i,0) = Y2(i,0) - A1(i);
// For each observation after the first...
for(j=1; j<Ny[1]; j++) {
  // For each time series...
  for(i=0; i<Ny[0]; i++) {
    // The current observation minus the drift...
    sum = Y2(i,j) - A1(i);
    // ...minus the dot product for the VAR(1) term...
    for(k=0; k<Ny[0]; k++)
      sum -= B2(i,k) * Y2(k,j-1);
    // ...is the innovation that yielded this observation
    V2(i,j) = sum;
  }
}



// VAR.forward
// 
// Supplied variables
// ----------------------------------------------------------------------------
// y    Modeled values (initially empty), [pn] array
// v    Innovation values, [pn] array
// a    Drift, [p] vector
// b    Lag-1 cross-correlation [pp] array

int i, j, k;
double sum;

// The first modeled value for each time series is just the first innovation
// plus the drift
for(i=0; i<Nv[0]; i++)
  Y2(i,0) = V2(i,0) + A1(i);
// For each innovation after the first...
for(j=1; j<Nv[1]; j++) {
  // For each time series...
  for(i=0; i<Nv[0]; i++) {
    // The drift plus the current innovation...
    sum = A1(i) + V2(i,j);
    // ...plus the dot product for the VAR(1) term...
    for(k=0; k<Nv[0]; k++)
      sum += B2(i,k) * Y2(k,j-1);
    // ...is the modeled output
    Y2(i,j) = sum;
  }
}



// Model.support
// ----------------------------------------------------------------------------


// Model parameters
struct parameters
{
  //--- Model parameters, direct ---
  PyArrayObject* d;  // Multivariate log-volatility offset
  PyArrayObject* e;  // VAR(1) coefficient, volatility shock to log-volatility
  PyArrayObject* g;  // Correlation, multivariate normal z vs. multivariate
                     // normal for innovation shocks

  //--- Model parameters, derivations ---
  PyArrayObject* rz; // Cross-correlation matrix, multivariate normal z to
                     // volatility shock v
  PyArrayObject* ri; // Inverted cross-correlation matrix, innovation shocks
  PyArrayObject* sc; // Scaling constants for multivariate normal pdf, given g
};


// An evaluation sample
struct sample
{
  // Log-likelihood of the time-series sample in the current evaluation, given
  // h at this point
  double Lxk;
  // Log-likelihood of the innovations up to the current one, given the history
  // of h to this point
  double Lx;
  // Multivariate normal proposal, which is just z for the current evaluation,
  // and for all time-series samples after the first in each log-volatility
  // computation
  double *z;
  // Multivariate log-volatility value based on the current evaluation, or the
  // previous proposal to use this struct if the current proposal hasn't been
  // crunched yet
  double *h;
  // Scratchpad 1-D arrays of the same length as z, h
  double *x0, *x1;
};


// Fast approximate exponent function EXP
//
// From N. Schraudolph, "A Fast, Compact Approximation of the Exponential
// Function," Neural Computation 11, 853-862 (1999)
static union {
  double d;
  struct {
#ifdef LITTLE_ENDIAN
    int j, i;
#else
    int i, j;
#endif
  } n;
} _eco;
#define EXP_A (1048576/M_LN2)
#define EXP_C 60801

#define EXP(y) (_eco.n.i = EXP_A*(y) + (1072693248 - EXP_C), _eco.d)
 

// Uniformly distributed random variate in [a,b)
double uniform(double a, double b)
{
  // Seed the PRNG on the first run only
  static int seeded = 0;
  if (!seeded) {
    int k;
    unsigned short seed[3];
    long int timeval = time(0);
    for(k=0; k<3; k++) {
      seed[k] = (unsigned short) timeval;
      timeval = timeval >> 16;
    }
    seed48(seed);
    seeded = 1;
  }
  // Compute and return a variate from U(a,b)
  return (b - a) * drand48() + a;
}


// Log probability density of the unit normal distribution
double normLogDensity(double x)
{
  static double logC = -log(2 * M_PI);
  return -0.5*x*x + logC;
}


// Random variate from a truncated normal distribution
double truncnorm(double a, double b)
{
  int k;
  register double x;
  double c, cp;
  const double diff = b-a;
  // Compute the scale value c as the maximum possible pdf (regular normal
  // distribution) within the truncated range
  if(a < 0) {
    if(b > 0) {
      // Center is within the range, so use the maximum pdf at center
      c = 1.02; // Higher than 1.0 due to max EXP +error of 1.02
    } else {
      // Range is below center, so use maximum pdf at its right edge
      c = 1.0625 * EXP(-0.5*b*b); // Scaled up by 1.02 / 0.96
    }
  } else {
    // Range is above center, so use pdf at its left edge
    c = 1.0625 * EXP(-0.5*a*a); // Scaled up by 1.02 / 0.96
  }
  // Accept-reject sampling for the truncated normal distribution
  do {
    // Generate a proposal within the truncated range
    x = a + diff*drand48();
    // Compute the pdf of the proposal, using the regular normal distribution
    cp = EXP(-0.5*x*x);
    // Accept (and exit the loop) when the scaled uniform variate is less than
    // the proposal's pdf
  } while(c*drand48() > cp);
  return x;
}


// Matrix multiplication x*y -> z, where the square matrix 'x' is a
// PyArrayObject and 'y' and 'z' are 1-D arrays.
void matrixMultiply(PyArrayObject* x, double *y, double *z)
{
  int kr, km;
  register double sum;
  const int rows = (int) x->dimensions[0];
  // Compute the dot product for each row of x
  for(kr=0; kr<rows; kr++) {
    // For this row of y and z...
    sum = 0;
    // Matrix multiplication dot product
    for(km=0; km<rows; km++)
      sum += PP2(x, kr, km) * PA1(y, km);
    PA1(z, kr) = sum;
  }
}


// Computes the log-volatility, given a proposal on z, the previously computed
// log-volatility value previous to the sample, and the model parameters.
void logVol(struct sample *S, struct parameters *P)
{
  int k;
  const int N = (int) P->d->dimensions[0];

  // Run the VAR(1) process to get the current multivariate log-volatility
  // value

  // Start with the matrix product for the VAR(1) term...
  matrixMultiply(P->e, S->h, S->x0);
  // ...then compute the multivariate volatility shock given the proposed
  // normal variate...
  matrixMultiply(P->rz, S->z, S->x1);
  // ...then add them, plus the multivariate volatility offset, to be the new h
  for(k=0; k<N; k++)
    SPA1(S, h, k) = SPA1(S, x0, k) + SPA1(S, x1, k) + SPP1(P, d, k);
}


// Computes the log-likelihood of a multivariate innovation x given a 1-D array
// containing a multivariate log-volatility.
void logLikelihood(struct sample *S, struct parameters *P, double *x)
{
  register int k;
  double val;
  const int N = (int) P->g->dimensions[0];
  
  // Inverse-scale the innovation by the square root of the volatilities in
  // preparation for decorrelation

  for(k=0; k<N; k++)
    // We could use a stripped-down version of exp here for speed; doing it
    // with the stdlib exp function takes about 70 CPU cycles (20 nS) on an
    // Intel QX9760 running at 3.8 GHz, though some of that is the unavoidable
    // array lookups and products in the equation.
    //
    SPA1(S, x0, k) = exp(-0.5 * SPA1(S, h, k)) * PA1(x, k);
    // Using EXP is very fast, but has error of about +/-6%. In a multivariate
    // test, the resulting error was -1112.29 vs an expected -1135.15.
    //SPA1(S, x0, k) = EXP(-0.5 * SPA1(S, h, k)) * PA1(x, k);

  // Now decorrelate the equi-variance innovations in preparation for the
  // log-likelihood computation for this multivariate log-volatility proposal
  matrixMultiply(P->ri, S->x0, S->x1);
  
  // Now compute log Pr(x[t]|h[t]), the log-likelihood of the decorrelated
  // multivariate innovation for this sample, conditional on the sample's
  // log-volatility.

  S->Lxk = 0;

  // For each time series...
  for(k=0; k<N; k++) {
    // Offset the decorrelated innovation by the mean, which is the correlation
    // between the normal variates for innovation and proposed volatility
    // shocks, times the value of the normal variate value underlying the
    // volatility shock for this sample
    SPA1(S, x0, k) = SPA1(S, x1, k) - SPP1(P, g, k) * SPA1(S, z, k);
    // The log-likelihood for this sample is simply the argument of the
    // exponential of the normal pdf, with the reciprocal scaling of the
    // exponential being done in log space by subtraction of the log of the
    // scaling coefficient
    val = SPP2(P, sc, k, 0) * SPA1(S, x0, k);
    // Fast squaring instead of pow()
    S->Lxk -= 0.5 * val*val + SPP2(P, sc, k, 1);
    // Since the innovation has been inverse-scaled by the square root of the
    // proposed volatility, the probability density needs to be scaled as
    // well. In log space, this is a subtraction of half the log-volatility.
    S->Lxk -= 0.5 * SPA1(S, h, k);
  }
  // Add the current log-likelihood to the accumulated log-likelihoods
  S->Lx += S->Lxk;
}



// Model.hybridGibbs
// 
// Supplied variables
// ----------------------------------------------------------------------------

//--- Supplied variables ---
// z    Independent normal variates, [pn] array (updated)
// x    Innovation values, [pn] array
// w    Wiggle value for +/- proposals (float)
// ni   Number of hybrid-Gibbs iterations (int)

//--- Model parameters, direct ---
// d    Volatility offset, [p] vector
// e    Lag-1 cross-correlations for VAR of volatilities, [pp] array
// g    Volatility/innovation shock correlations, [p] vector

//--- Model parameters, derivations ---
// rz   Concurrent xcorr, innovation-volatility normals, [pp]
// ri   Inverse concurrent xcorr, innovation shocks, [pp]
// sc   Constants for multivariate normal PDF


// Run evaluations until odds of error in selection between current and
// proposed z will be less than 1/100.
#define MIN_DIFF 0.01

// The number of time series making up a single multivariate sample
const int Nt = Nx[0];
// The number of multivariate samples
const int Ns = Nx[1];

int ki, ke;
int k0, k1, k2, k3;
double val;

// Multivariate innovation for one current time-series sample
double *xk;
// The multivariate proposal on z
double zp[Nt];
// The multivariate value of h for current and proposal evaluations, at the
// first time-series sample of the evaluation interval
double he[2][Nt];

// A set of evaluation sample structs
struct sample se[2];
// A struct for the model parameters
struct parameters P;


// Generate parameter struct
P.d = d_array; P.e = e_array; P.g = g_array;
P.rz = rz_array; P.ri = ri_array; P.sc = sc_array;

// Initialize various arrays
xk = (double *) malloc(sizeof(double) * Nt);
for(k0=0; k0<2; k0++) {
  se[k0].z = (double *) malloc(sizeof(double) * Nt);
  se[k0].h = (double *) malloc(sizeof(double) * Nt);
  se[k0].x0 = (double *) malloc(sizeof(double) * Nt);
  se[k0].x1 = (double *) malloc(sizeof(double) * Nt);
}

// For each iteration...
for(ki=0; ki<ni; ki++) {

  // The "previous" log-volatility of the first time-series sample is set to
  // the log-volatility offset plus a simulated, latent-parameter value.
  // TODO
  for(k0=0; k0<2; k0++) {
    for(k1=0; k1<Nt; k1++)
      PA1(se[k0].h, k1) = D1(k1);
  }

  // For each time-series sample...
  for(k0=0; k0<Ns; k0++) {
    
    // Simulate two sets of h samples from the current z and a proposal on z,
    // with a likelihood evaluation that starts from the current time-series
    // sample to some subsequent time-series sample at which the log-likelihood
    // of the innovation for that sample becomes substantially the same for the
    // current and proposal value of z.
  
    k1 = k0;
  
    // Clear the log-likelihood and normal log-density accumulators of each
    // evaluation struct for the next pair of evaluations.
    for(k2=0; k2<2; k2++)
      se[k2].Lx = 0;
    
    do {
  
      // Prepare current-innovation vector xk
      for(k2=0; k2<Nt; k2++)
        PA1(xk, k2) = X2(k2, k1);
      
      // For each evaluation...
      for(k2=0; k2<2; k2++) {
  
        // Set z...
        for(k3=0; k3<Nt; k3++) {
          val = Z2(k3, k1);
          if(k1==k0) {
            if(k2==1) {
              // First time-series sample for the second evaluation, use a
              // proposal on z instead of z itself
              val = truncnorm(val-w, val+w);
              zp[k3] = val;
            }
          }
          PA1(se[k2].z, k3) = val;
        }
        // ...compute the multivariate log-volatility
        logVol(&se[k2], &P);
        // ...and the log-likelihood of the innovation for this time-series
        // sample, given the log-volatility
        logLikelihood(&se[k2], &P, xk);
        
        // Save the first log-volatility
        if(k1==k0) {
          for(k3=0; k3<Nt; k3++)
            he[k2][k3] = PA1(se[k2].h, k3);
        }
      }
  
      // Ready for the next time series sample, assuming there is one
      k1++;

      // Continue only if difference in log-likelihood of the current
      // time-series sample of the evaluations between current and proposed z
      // is great enough to warrant doing so.
      val = se[0].Lxk - se[1].Lxk;
    
    } while (k1<Ns && (val > MIN_DIFF || val < -MIN_DIFF));

    // Metropolis-hastings selection of current or proposal
    val = se[1].Lx  -  se[0].Lx;
    if(val > 0 || log(uniform(0, 1)) < val) {
      // Proposal accepted
      ke = 1;
      for(k2=0; k2<Nt; k2++)
        Z2(k2, k0) = zp[k2];
    } else {
      // Proposal rejected
      ke = 0;
    }
    
    // Set the initial log-volatility of the evaluations for the next
    // time-series sample to the log-volatility of the selected evaluation for
    // this time-series sample.
    for(k2=0; k2<Nt; k2++) {
      val = he[ke][k2];
      for(k3=0; k3<2; k3++)
        PA1(se[k3].h, k2) = val;
    }
  
    // Done with log-volatility draw for this time-series sample
  }

  // Done with iterations
}

// Free array memory
free(xk);
for(k0=0; k0<2; k0++) {
    free(se[k0].z);
    free(se[k0].h);
    free(se[k0].x0);
    free(se[k0].x1);
}



// Model.likelihood
// 
// Supplied variables
// ----------------------------------------------------------------------------

//--- Supplied variables ---
// z    Independent normal variates, m previously simulated sets, [mpn] array
// x    Innovation values, [pn] array

//--- Model parameters, direct ---
// d    Volatility offset, [p] vector
// e    Lag-1 cross-correlations for VAR of volatilities, [pp] array
// g    Volatility/innovation shock correlations, [p] vector

//--- Model parameters, derivations ---
// rz   Concurrent xcorr, innovation-volatility normals, [pp]
// ri   Inverse concurrent xcorr, innovation shocks, [pp]
// sc   Constants for multivariate normal PDF

//--- Return value (double float) ---
// The linear likelihood P(x|w) of the innovations given the parameters,
// integrated over the log-volatility simulations along the first axis of z


// The number of time series making up a single multivariate sample
const int Nt = Nx[0];
// The number of multivariate samples
const int Ns = Nx[1];

int ki;
int k0, k1;
double val, maxVal;

// Multivariate innovation for one current time-series sample
double *xk;
// Log-likelihood of the innovation for the time-series sample in the current
// evalution, accumulated for each iteration
double Lx[Nz[0]];

// A single evaluation sample struct
struct sample se;
// A struct for the model parameters
struct parameters P;


// Generate parameter struct
P.d = d_array; P.e = e_array; P.g = g_array;
P.rz = rz_array; P.ri = ri_array; P.sc = sc_array;

// Initialize various arrays
xk = (double *) malloc(sizeof(double) * Nt);
se.z = (double *) malloc(sizeof(double) * Nt);
se.h = (double *) malloc(sizeof(double) * Nt);
se.x0 = (double *) malloc(sizeof(double) * Nt);
se.x1 = (double *) malloc(sizeof(double) * Nt);


// Iterate over each set of simulated z values...
for(ki=0; ki<Nz[0]; ki++) {

  // Initialize this iteration's accumulator
  Lx[ki] = 0;

  // The "previous" log-volatility of the first time-series sample is set to
  // the log-volatility offset plus a simulated, latent-parameter value.
  // TODO
  for(k1=0; k1<Nt; k1++)
    PA1(se.h, k1) = D1(k1);

  // For each time-series sample...
  for(k0=0; k0<Ns; k0++) {
    
    // For each time series...
    for(k1=0; k1<Nt; k1++) {
      // Prepare current-innovation vector xk...
      PA1(xk, k1) = X2(k1, k0);
      // ...and set z
      PA1(se.z, k1) = Z3(ki, k1, k0);
    }
    
    // Compute the multivariate log-volatility...
    logVol(&se, &P);
    // ...and the log-likelihood of the innovation for this time-series
    // sample, given the log-volatility
    logLikelihood(&se, &P, xk);
    // Accumulate the log-likelihood
    Lx[ki] += se.Lxk;
  
    // Done with log-likelihood evaluation for this time-series sample
  }

  // Done with iterations
}

// Free array memory
free(xk);
free(se.z);
free(se.h);
free(se.x0);
free(se.x1);


// Compute the result, the integrated P(x|w) over all the sets of simulated z
// values

// Compute the maximum log density and save it to subtract from everybody so
// that the maximum log density is normalized to zero
maxVal = -HUGE_VAL;
for(ki=0; ki<Nz[0]; ki++) {
  if(Lx[ki] > maxVal)
    maxVal = Lx[ki];
}
// Accumulate the linearized normalized densities
val = 0;
for(ki=0; ki<Nz[0]; ki++)
  val += exp(Lx[ki] - maxVal);
// ...and return the denormalized, log mean of the linear densities...
return_val = log(val) + maxVal - log(Nz[0]);