liborprocess.java

金融资产定价,随机过程,MONTE CARLO 模拟 JAVA 程序和文档资料

/* WARANTY NOTICE AND COPYRIGHT
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
GNU General Public License 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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

Copyright (C) Michael J. Meyer

[e-mail]
[e-mail]

*/


/*
 * LiborProcess.java
 *
 * Created on August 29, 2002, 8:41 AM
 */

package Libor.LiborProcess;

import ArrayClasses.*;
import Statistics.*;
import LinAlg.*;
import Triggers.*;

/**<p>The basic class simulating Libor paths as well as the paths of two
 * log-normal approximations (<code>X0,X1</code>, see <i>LiborProcess.ps</i>).
 * Time steps move directly from one point <code>T_i</code> on the Libor
 * tenor structure to the following point <code>T_{i+1}</code> using a predictor 
 * corrector simulation. Paths are simulated in the forward martingale 
 * measure <code>P_n</code> at the terminal time <code>T=T_n</code>.</p>
 *
 * <p>This setup can handle all Libor derivatives with payoffs which are 
 * deterministic functions of the matrix <code>(L_j(T_k))_{j&lt;=k})</code>
 * of Libors evaluated at tenor points <code>T_k</code>. This covers caps,
 * swaptions, certain Bermudan swaptions etc. but not derivatives which 
 * depend on Libors evaluated more frequently.</p>
 *
 * <p>The class provides methods to compute Libor sample paths and basic 
 * functionals of the Libor path such as zero coupon bonds, annuity numeraires,
 * swap rates etc.</p>
 * 
 * <p>It also provides methods to compute log-normal approximations
 * to certain vectors of Libors which can be simulated directly and used for
 * fast approximate valuation of some Libor derivatives.</p>
 * 
 *
 * @author  Michael J. Meyer
 */
public class LiborProcess {
    
    int n;               // number of forward Libors
    double[] delta,      // delta[t]=T_{t+1}-T_t,length of compounding intervals
             Tc,         // tenor structure, Tc[t]=T_t
             l,          // initial Libors, l[j]=L_j(0)
             x;          // initial X-Libors x[j]=X_j(0)=delta_jL_j(0)
    
    // volatility and correlation structure
    FactorLoading factorLoading;
    
    // The array of Wiener increments (independent standard normal
    // increments) driving the discretized Libor paths.
    // Allocated as upper triangular array.
    double[][] Z;
    
    // The following are allocated as lower triangular arrays:
    
    // Array containing the X-Libor paths X_j(t)=delta_jL_j(t), for
    // t=T_0,T_1,...,T_j (discrete t=0,1,...,j). 
    // Allocating this as a lower triangular array allows
    // natural indices X[j][t]=X_j(t), t=T_k, k<=j.
    double[][] X;
    
    // Array containing the Y-Libor paths Y_j(t)=log(X_j(t)), for
    // t=T_0,T_1,...,T_j.
    double[][] Y;
    
    // Array containing the log-Gaussian approximation X^0_j(t)
    // of X_j(t). These paths are computed along with the Libor paths
    double[][] X0;
    
    // Array containing the Gaussian log Y0_j(t)=log(X^0_j(t))
    double[][] Y0;
    
    // Array containing the log-Gaussian approximation X^1_j(t)
    // of X_j(t). These paths are computed along with the Libor paths
    double[][] X1;
    
    // Array containing the Gaussian log Y0_j(t)=log(X^0_j(t))
    double[][] Y1;
    
    double[] alpha, beta;    // linearization constants for X^1

    // drift step vector
    double[] m;
    
    // volatility step vector
    double[] V;
    
    // vector used to store factors X_j(t)/(1+X_j(t)) during time steps
    double[] F;
    
    // Reference to the data array containing the sequence 
    // {@link FactorLoading#covariationMatrices}. Direct access to the data 
    // array to gain speed.
    double[][][] covariationMatrices;
    
    // Reference to the data array containing the sequence 
    // {@link FactorLoading#choleskyRoots} of Cholesky roots of the
    // covariation matrices. Direct access to the data array to gain speed.
    double[][][] choleskyRoots;
    

/*******************************************************************************
 *
 *                 ACCESSORS, FORWARDING TO FACTORLOADING
 *
 ******************************************************************************/ 
    
    /** The number <code>n</code> of Libors including <code>L_0</code>.
     */
    public int dimension(){ return n; }
    
    /** The array <code>Tc[j]=T_j</code> (continuous times).
     */
    public double[] tenorStructure(){ return Tc; }
    
    
    /** The array <code>delta[j]=delta_j</code> of accrual periods.
     */
    public double[] delta(){ return delta; }
    
    
    /** <p>The deterministic volatility <code>sigma_i(t)</code> of Libor
     *  <code>L_i(t)</code>.</p>
     *
     * @param i Libor index.
     * @param t current discrete time.
     */
    public double sigma(int i, int t)
    { return factorLoading.sigma(i,t); }

    
    /** <p>The instantaneous correlations of Libor increments
     *  <code>dL_i,dL_j</code>. Recall that these are assumed to be constant.
     *  </p>
     *
     * @param i Libor index.
     * @param j Libor index.
     */
    public double rho(int i, int j){ return factorLoading.rho(i,j); } 

    
    /** The array <code>l[j]=L_j(0)</code> of initial Libors.
     */
    public double[] initialTermStructure(){ return l; }
    
    
    /** The array <code>x[j]=X_j(0)</code> of initial X-Libors.
     */
    public double[] xInitial(){ return x; }
    
    
    /** The {@link FactorLoading} of the Libor process
     */
    public FactorLoading factorLoading(){ return factorLoading; }
    
    
   /** <p>The <code>n by n</code> matrix of instantaneous log-Libor 
    *  correlations <code>(rho_ij)_{0&lt;=i,j&lt;n}</code>. Note that this 
    *  includes Libor <code>L_0</code> ie. we are taking the view here that
    *  all Libors live to the horizon. See document <i>LiborProcess.ps</i>.</p>
    */
    public ColtMatrix correlationMatrix()
    {
        return factorLoading.rho();
    }
    
    
    /** Prints the matrix <code>(rho_ij)</code> of instantaneous 
     *  log-Libor correlations for 
     */
    public void printCorrelationMatrix()
    {
        System.out.println(factorLoading.rho().toString());
    }
    
    
  /** <p>The log-Libor covariation-matrix
   * <p>
   * <center>
   * <code>CV(p,q,t,T)=( integral_t^T cv_ij(s)ds )_{i,j=p}^{q-1},</code><br
   * where<br>
   * <code>cv_ij(s)=sigma_i(s)sigma_j(s)rho_ij=nu_i(s).nu_j(s).</code>
   * </center>
   * </p>
   *
   * <p>WARNING: index shift to zero based matrix indices 
   * <code>(i,j)-&gt(i-p,j-p)</code>.</p>
   *
   * @param p Libors <code>L_j, j=p,...,q-1.</code>
   * @param q Libors <code>L_j, j=p,...,q-1.</code>
   * @param t time interval <code>[t,T]</code>, continuous time.
   * @param T time interval <code>[t,T]</code>, continuous time.
   */ 
   public ColtMatrix 
   logCovariationMatrix(int p,int q, double t, double T)
   {
       return factorLoading.logCovariationMatrix(p,q,t,T);
   }
   

  /** <p>The Cholesky root of 
   *  {@link #logCovariationMatrix(int,int,double,double)}.</p>
   *
   * @param p Libors <code>L_j, j=p,...,q-1.</code>
   * @param q Libors <code>L_j, j=p,...,q-1.</code>
   * @param t time interval <code>[t,T]</code>, continuous time.
   * @param T time interval <code>[t,T]</code>, continuous time.
   */ 
   public ColtMatrix 
   logCovariationCholeskyRoot(int p,int q, double t, double T)
   {
       return factorLoading.logCovariationCholeskyRoot(p,q,t,T);
   }
    
/*******************************************************************************
 *
 *                             CONSTRUCTOR
 *
 ******************************************************************************/    

    /** Function <code>f(y)=exp(y)/(1+exp(y))</code> which is linearized
     *  for the <code>X1</code>-dynamics and needed to define the linearization
     *  constants <code>alpha_k, beta_k</code>.
     */
    private double f(double y)
    {
        double x=Math.exp(y);
        return x/(1+x);
    }
    
    
    /**
     * @param l array of initial Libors <code>l[i]=L_i(0)</code>.
     * @param fl volatility and correlation structure, see
     * {@link FactorLoading}.
     */
    public LiborProcess(double[] l, FactorLoading fl) {
        
        n=fl.getDimension();
        Tc=fl.getTenorStructure();
        this.l=l;                       // initial term stucture
        factorLoading=fl;
        
        // allocate delta array and compute deltas
        delta=new double[n];
        for(int j=0;j<n;j++)delta[j]=Tc[j+1]-Tc[j];
        
        // initial X-Libors X_j(0)=delta_jL_j(0)
        x=new double[n];
        for(int j=0;j<n;j++)x[j]=delta[j]*l[j];

        
        // allocate array for standard normal increments,
        // direct access to data for more speed.
        Z=new UpperTriangularArray(n-1).getData();
        
        // allocate path arrays, direct access to data for more speed:
        
        // path array, X-dynamics, Y_j=log(X_j)
        X=new LowerTriangularArray(n).getData();
        Y=new LowerTriangularArray(n).getData();
        
        // path array, X0-dynamics
        X0=new LowerTriangularArray(n).getData();
        Y0=new LowerTriangularArray(n).getData();
        
        // path array, X1-dynamics
        X1=new LowerTriangularArray(n).getData();
        Y1=new LowerTriangularArray(n).getData();
        
        // initialize path arrays
        for(int j=0;j<n;j++)
        { 
            X[j][0]=x[j];      Y[j][0]=Math.log(X[j][0]);
            X0[j][0]=X[j][0];  Y0[j][0]=Math.log(X0[j][0]); 
            X1[j][0]=X[j][0];  Y1[j][0]=Math.log(X1[j][0]);
        }
            
        m=new double[n];
        V=new double[n];
        F=new double[n];
        
        covariationMatrices=
        factorLoading.getCovariationMatrixArray().getMatrices();
        
        choleskyRoots=
        factorLoading.getCholeskyRootArray().getMatrices();
        
        // allocate and initialize linearization constants alpha_k, beta_k
        // for the X1-dynamics (includes X1_0, no index downshift j->j-1)
        alpha=new double[n]; beta=new double[n];
        for(int j=0;j<n;j++){
            
            double a=Y[j][0]-0.6*Math.log(3),
                   b=Y[j][0]+0.4*Math.log(3);
            
            alpha[j]=(b*f(a)-a*f(b))/(b-a);
            beta[j] =(f(b)-f(a))/(b-a);
        }
        
    } // end constructor
    
    
    /** Constructs Libor process from {@link LMM_Parameters} parameter
     *  object.
     *
     * @param lmmParams parameter object.
     */
    public LiborProcess(LMM_Parameters lmmParams)
    {
        this(lmmParams.initialTermStructure(),   
             lmmParams.factorLoading());
    }

    
/*******************************************************************************
 *
 *                  WIENER INCREMENT GENERATION
 *
 ******************************************************************************/ 
    
    
    /** <p>Fills the <code>Z-array</code> with a new set of independent standard
     *  normal increments needed to evolve Libors from discrete time 
     *  <code>t</code> to discrete time <code>T</code>. It does not matter if 
     *  the whole set of Libors or a subset is evolved, the same set of 
     *  Z-increments has to be generated.
     *
     * @param t begining of simulation.
     * @param T end of simulation.
     */
    public void newWienerIncrements(int t, int T)
    {
        for(int k=t;k<T;k++)
        for(int j=k+1;j<n;j++)Z[k][j-k-1]=Random.STN();
     
    } // end newWienerIncrements
        

/*******************************************************************************
 *
 *                         TIME STEP
 *
 ******************************************************************************/    

    /* Remember that we compute paths with time steps T_t->T_{t+1}, that is, we
     * are stepping from one Libor reset time to the next. Discrete (integer)
     * time <code>t</code> corresponds to continuous time 
     * <code>Tc[t]=T_t</code>.
     *
     * Recall also that X_j(t)=delta_jL_j(t) and Y_j(t)=log(X_j(t)).
     * With the above conventions about discrete time we have X_j(T_j)=X[j][j].
     * 
     * We simulate the path of the driftless approximation along with the Libor 
     * path since this allows us to get highly correlated control variates but 
     * adds only a very small computational overhead.
     */
    

     /** <p>Evolves the X-Libors <code>L_j(t), j&gt=p</code> from time 
      *  <code>T_t</code> (discrete time <code>t</code>) to time 
      *  <code>T_{t+1}</code> (discrete time <code>t+1</code>) in a single 
      *  time step.</p>
      *
      * <p> Assumes that the array Z of driving Wiener increments has been
      * filled with new values. See document 
      * <i>LiborProcessImplementation.ps</i></p>
      *
      * <p>Any of the three dynamics <code>X,X0,X1</code> (see document
      * <i>Liborprocess.ps</i>, section Implementation notes) can be moved
      * forward. Set the corresponding boolean flags to indicate which.</p>
      * 
      * @param t current discrete time (continuous time <code>T_t</code>).
      * @param p Libors evolved are <code>L_j, j=p,p+1,...,n-1.</code>
      * @param Xpath do the time step for the X-dynamics (yes/no).
      * @param X0path do the time step for the X0-dynamics (yes/no).
      * @param X1path do the time step for the X1-dynamics (yes/no).
      */
     public void timeStep
     (int t, int p, boolean Xpath, boolean X0path, boolean X1path)
     {
         /* Ordinarily a method such as this would be broken up into smaller
          * methods (handling X,X0,X1 separately). However this method will
          * be called in loops and so we want to save the overhead of calling 
          * other functions.
          */
         
         double[][] C=covariationMatrices[t], 
                    L=choleskyRoots[t];
                    
         int q=Math.max(t+1,p);   
         // only Libors L_j, j>=q make the step. To check the index shifts 
         // to zero based java array indices below consider the case p=t+1 
         // (all Libors), then q=t+1.
         
         
         // volatility step vector V, common to all dynamics X0,X1,X
         for(int i=q;i<n;i++)
         {
             V[i]=0;
             // note shift to zero based index in L[t] is j->j-t-1.
             for(int j=t+1;j<=i;j++) V[i]+=L[i-t-1][j-t-1]*Z[t][j-t-1];  
         }
    
         
         // X0-DYNAMICS
         
         if(X0path){
           // drift step vector for Y^0(t)=log(X^0(t)).
           // C is upper triangular
           // shift to zero based indices for C: (i,j)->(i-t-1,j-i)
           // so main diagonal has j-i=0.
           for(int j=q+1;j<n;j++) F[j]=1-1/(1+x[j]);
           for(int i=q;i<n;i++){ 
             
               m[i]=-0.5*C[i-t-1][0];
               for(int j=i+1;j<n;j++)m[i]-=F[j]*C[i-t-1][j-i]; 
           } // end for i
         
           // evolve the log of the driftless path
           for(int i=q;i<n;i++){ Y0[i][t+1]=Y0[i][t]+m[i]+V[i];