OrnsteinUhlenbeckWithIntegratedProcess.java

package umontreal.ssj.stochprocess;
 
import umontreal.ssj.randvar.*;
import umontreal.ssj.randvarmulti.*;
// import umontreal.ssj.stochprocess.OrnsteinUhlenbeckProcess;
import cern.colt.matrix.impl.*;

public class OrnsteinUhlenbeckWithIntegratedProcess extends OrnsteinUhlenbeckProcess {
   double[] integratedPath; // length [d + 1]
   MultinormalGen[] normalsCorrGen; // length [d]
   double[] bdt; // length [d]: (beta * dt)
   double[] oneMExpMADtOverAlpha; // length [d]: (1.0 - exp(-alpha * dt)) / alpha
   double c22; // sigma * sigma / 2.0 / alpha^3

   public OrnsteinUhlenbeckWithIntegratedProcess(double x0, double alpha, double b, double sigma, NormalGen gen) {
      super(x0, alpha, b, sigma, gen);
   }

   @Override
   protected void initArrays(int d) {
      super.initArrays(d); // [Not sure why the "d" is there, but was there in super.initArrays...]
      // this is called by init() which is called in
      // StochasticProcess.setObservationTimes().
      integratedPath = new double[d + 1];
      bdt = new double[d];
      oneMExpMADtOverAlpha = new double[d];
      // iTime == d not included
      for (int iTime = 0; iTime < d; iTime++) {
         double dt = t[iTime + 1] - t[iTime];
         bdt[iTime] = beta * dt;
         double oneMExpMADt = -Math.expm1(-alpha * dt); // (1.0 - exp(-alpha * dt);
         oneMExpMADtOverAlpha[iTime] = oneMExpMADt / alpha;
      }
 
      double c11 = sigma * sigma / 2.0 / alpha;
      double c12 = c11 / alpha;
      c22 = c12 / alpha;
      normalsCorrGen = new MultinormalGen[d];
      for (int iTime = 0; iTime < d; iTime++) {
         double dt = t[iTime + 1] - t[iTime];
         double expMADt = alphadt[iTime];
         double oneMExpMADt = -Math.expm1(-alpha * dt); // (1.0 - exp(-alpha * dt);
 
         DenseDoubleMatrix2D covar = new DenseDoubleMatrix2D(2, 2);
         covar.set(0, 0, c11 * oneMExpMADt * (1.0 + expMADt)); // sigma^2/2/alpha * (1 - exp(-2 * alpha * dt))
         covar.set(0, 1, c12 * oneMExpMADt * oneMExpMADt); // sigma^2/2/alpha^2 * (1 - exp(-alpha*dt))^2
         covar.set(1, 0, covar.get(0, 1));
         covar.set(1, 1, c22 * (-3.0 + 2.0 * alpha * dt + 4.0 * expMADt - expMADt * expMADt));
 
         double[] mu = { 0.0, 0.0 }; // the expectation terms depend on the previous path values and are added
                                     // elsewhere.
         normalsCorrGen[iTime] = new MultinormalPCAGen(gen, mu, covar);
      }
   }

   @Override
   public double[] generatePath() {
      path[0] = x0;
      integratedPath[0] = 0.0;
      double[] normalsCorr = new double[2];
      for (int iTime = 1; iTime <= d; iTime++) {
         normalsCorrGen[iTime - 1].nextPoint(normalsCorr);
         path[iTime] = badt[iTime - 1] + path[iTime - 1] * alphadt[iTime - 1] + normalsCorr[0];
 
         integratedPath[iTime] = integratedPath[iTime - 1] + bdt[iTime - 1]
               + oneMExpMADtOverAlpha[iTime - 1] * (path[iTime - 1] - beta) + normalsCorr[1];
      }
      observationIndex = d;
      return path;
   }

   public double[] getIntegratedPath() {
      return integratedPath;
   }

   public double[] getExpectedFutureDiscount(double[] couponTimes) {
      double[] analyticNotificationDiscounts = new double[couponTimes.length];
      for (int iTime = 1; iTime < couponTimes.length; iTime++) {
         int iTimePath = (iTime < d) ? iTime : d;
         double dt = couponTimes[iTime] - t[iTimePath]; // noticePeriod
         double expMADt = Math.exp(-alpha * dt);
         double mu_2 = beta * dt + (1.0 - expMADt) / alpha * (path[iTimePath] - beta);
         double sigma22 = c22 * (-3.0 + 2.0 * alpha * dt + 4.0 * expMADt - expMADt * expMADt);
         analyticNotificationDiscounts[iTime] = Math.exp(-mu_2 + sigma22 / 2.0 - integratedPath[iTimePath]);
      }
      return analyticNotificationDiscounts;
   }

   public double[] getTotalAnalyticDiscount(double[] times) {
      double[] analyticDiscount = new double[times.length];
      analyticDiscount[0] = 1.0;
      for (int iTime = 1; iTime < times.length; iTime++) {
         double dt = times[iTime] - times[0];
         double expMADt = Math.exp(-alpha * dt);
         double mu_2 = beta * dt + (1.0 - expMADt) / alpha * (x0 - beta);
         double sigma22 = c22 * (-3.0 + 2.0 * alpha * dt + 4.0 * expMADt - expMADt * expMADt);
         analyticDiscount[iTime] = Math.exp(-mu_2 + sigma22 / 2.0);
      }
      return analyticDiscount;
   }
}