liborprocess.java

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

             
             public double getValue(int t)
             {
                 newPathSegment(t,j,j); return Math.log(L1(j,j));
             }
         }; // end return new
         
     } // end logL1

     
             
/*******************************************************************************
 *
 *                   LOG-GAUSSIAN APPROXIMATE LIBOR 
 *
 ******************************************************************************/
     
     /** <p>The vector <code>(X^0_p(T_p),...,X^0_{n-1}(T_p)</code>, that is,
      *  snapshot of the log-normal <code>X0</code> Libors 
      *  <code>X^0_j, j=p,p+1,...,n-1</code> at time <code>T_p</code>.</p>
      *  
      *  This vector is highly correlated to the true Libors at time 
      *  <code>T_p</code> and allows very fast direct simulation 
      *  (the logarithms are multinormal). 
      *  Used for log-normal approximate pricing and in control variates.</p>
      *
      * <p>WARNING: no conditioning, time <code>t</code> ignored,
      * sampling corresponds to Libor path computation from time 
      * <code>t=0</code>.</p>
      *
      * @param p Libors <code>X_j(T_p), j&gt;=p</code>.
      * Must satisfy <code>p&lt;n</code>.
      */
     public LiborVector X0LiborVector(final int p)
     {   
         final int N=n;     // avoid ambiguity 
         
         // if p=0 the vector is the constant X(0)
         if(p==0)
         return new LiborVector(this,p){
             
             // write next sample vector into the array U
             public void nextSample(){ for(int i=0;i<N;i++)U[i]=X(i,0); }

         }; // end return new
                 
         // if p>0
         
         final double[] Z=new double[n-p];  // standard normal vector
         final ColtMatrix covariationMatrix=new ColtMatrix(n-p,n-p);
         
         // set covariations
         for(int i=p;i<n;i++)
         for(int j=p;j<n;j++){
             
             double cv_ij=factorLoading.integral_sgi_sgj_rhoij(i,j,0,Tc[p]);
             covariationMatrix.setQuick(i-p,j-p,cv_ij);
             covariationMatrix.setQuick(j-p,i-p,cv_ij);
         }
         
         final ColtMatrix choleskyRoot=covariationMatrix.choleskyRoot();
         
         // use the raw data arrays for faster simulation
         final double[][] chr=choleskyRoot.toArray();
         //System.out.println(ArrayUtils.toString(chr));
         
         final double[] // means of the log(X^0_j(T_m))
                        logMeans=new double[n-p];
         
         // set means of logarithms
         for(int i=p;i<n;i++){
             
            logMeans[i-p]=Math.log(x[i])-
                          0.5*factorLoading.integral_sgi_sgj_rhoij(i,i,0,Tc[p]);
            double F=x[i]/(1+x[i]), 
                   sum=0;
            for(int j=i+1;j<n;j++)
            sum+=factorLoading.integral_sgi_sgj_rhoij(i,j,0,Tc[p]);
            sum*=F;
            logMeans[i-p]-=sum;
         }
         
         
         // the Libor vector
         LiborVector X0p=new LiborVector(this,p){
             
             // write next sample vector into the array U
             public void nextSample()
             {
                 for(int i=p;i<N;i++)Z[i-p]=Random.STN();
                 for(int i=p;i<N;i++)
                 {
                     // log-Libor log(X^0_i(T_m))
                     double y_i=logMeans[i-p];
                     for(int j=p;j<=i;j++)y_i+=chr[i-p][j-p]*Z[j-p];
                     
                     U[i-p]=Math.exp(y_i);
                 }
             } // end nextSample()
         }; // end X0p
         
         return X0p;
      
     } // driftlessLibors
                     
             
/*******************************************************************************
 *
 *                          ZERO COUPON BONDS 
 *
 ******************************************************************************/

        
 
    /** The zero coupon bond <code>B_i(0)=B(0,T_i)</code>.
     *  
     * @param i bond matures at time <code>T_i<code>
     */
     public double B0(int i)
     { 
         double f=1;
         // accumulate 1 from time t=0 to time t=T_i                    
         for(int j=0;j<i;j++)f*=(1+x[j]); 
         return 1/f;                                 // B_i(0)=1/f  
     }   
 
    /** The zero coupon bond <code>B_i(T_t)=B(s,T)</code> with 
     *  <code>s=T_t,T=T_i</code>. Assumes that all Libors 
     *  <code>L_k, k=t,t+1,...,n-1</code> have been computed up to time 
     *  <code>s=T_t</code> (discrete time <code>t</code>).
     *
     *  @param t bond evaluated at time <code>T_t</code>
     *  @param i bond matures at time <code>T_i</code>, 
     *  must satisfy <code>t&lt=i</code>.
     */
    public double B(int i, int t)
    {     
        double f=1; 
    
        // accumulate 1 forward from time t to time i
        for(int k=t;k<i;k++)f*=(1+X[k][t]); 
        return 1/f;                                 // B_i(T_t)=1/f  
    }

    /** <p>Accrual factor <code>B_i(t)/B_n(t)</code>. This factor shifts
     *  a cashflow from time <code>T_i</code> to the horizon <code>T_n</code>
     *  with Libors in state at time <code>T_t</code> (discrete time
     *  <code>t</code>). In other words the shift is carried out at discrete
     *  time <code>t</code>.
     *  
     *  <p>Needs all Libors <code>L_j, j&gt=i</code>
     *  up to time <code>T_t</code> (discrete time <code>t</code>).</p>
     *
     * @param t current discrete time (continuous time <code>T_t</code>).
     * @param i cashflow shifted from time <code>T_i</code> to horizon.
     */
     public double forwardTransport(int t, int i)
     { 
         double f=1.0;
         for(int j=i;j<n;j++)f*=(1+X[j][t]);
         return f;
     } 
     
     
    /** <p>Accrual factor <code>1/B(T_i,T_n)=1/B(n,i)</code>. This factor shifts
     *  a cashflow from time <code>T_i</code> to the horizon <code>T_n</code>
     *  with Libors in state at time <code>T_i</code>.</p>
     *  
     *  Needs all Libors <code>L_j, j&gt=i</code>
     *  up to time <code>T_i</code> (discrete time <code>i</code>).</p>
     *
     * @param i cashflow shifted from time <code>T_i</code> to horizon.
     */
     public double forwardTransport(int i)
     { 
         if(i==n)return 1.0;
         
         double f=1.0;
         for(int j=i;j<n;j++)f*=(1+X[j][i]);
         return f;
     } 
     
     /** Discrete time <code>j</code> such that 
      *  <code>T_j<&lt=t&lt;T_{j+1}</code>.
      *
      * @param t continuous time
      */
     private int j(double t)
     {
         int j=0, k=n, i;
         // binary search
         while(k-j>1){  i=(j+k)/2; 
                        if(t<Tc[i]) k=i; else j=i; }
         
         return j;
     } // end j
     

     /** <p>The zero coupon bond <code>B(t,T)</code> with <code>t&lt=T</code>
      *  arbitrary times. Simple discounting on subintervals of Libor
      *  accrual intervals.</p>
      *
      * <p>Note: since Libors are not available at time <code>t</code> they are
      *  evaluated at the time <code>T_i</code> closest to <code>t</code>.
      *  Consequently this is only an approximation to the zero coupon bond
      *  <code>B(t,T)</code>. In fact this bond can be computed accurately
      *  in the LMM only if both <code>t</code> and <code>T</code> coincide with
      *  Libor reset times.</p>
      *
      * <p>Needs all Libors <code>L_j, j&gt=i</code>
      * up to time <code>T_{i+1}</code> where <code>T_i&lt=t<T_{i+1}</code>.
      * </p>
      *
      * @param t current continuous time.
      * @param T continuous time of bond maturity, must satsisfy
      * <code>t&lt=T</code>.
      */
     public double B(double t, double T)
     {
         int i=j(t),    // T_i<=t<T_{i+1}
             j=j(T);    // T_j<=T<T_{j+1}
         
         // time T_u closest to t
         int u=(t-Tc[i])<(Tc[i+1]-t)? i : i+1;
         
         double f;              // forward accrual factor
         
         if(i==j)  // case i=j, ie. T_i<=t<=T<T_{i+1}
         {
             // Libor already set at time T_i
             f=1+(T-t)*X[i][i]/delta[i];
         }
         else      // case i<j, ie. T_i<=t<T_{i+1}<=T_j<=T
         {
             // accrual factor from time t to time T_{i+1}, 
             // Libor already set at time T_i
             f=1+(Tc[i+1]-t)*X[i][i]/delta[i];
             // acrual forward from time T_{i+1} to time T_j
             // Libor evaluated at time T_u as a proxy for t
             for(int k=i+1;k<j;k++)f*=(1+X[k][u]);
             // acrual forward from time T_j to time T
             f*=1+(T-Tc[j])*X[j][u]/delta[j];
         }
         
         return 1.0/f;
         
     } // end B(t,T)

     
/*******************************************************************************
 *
 *                        FORWARD SWAP RATES
 *
 ******************************************************************************/                         
             

   /** <p>The forward swap rate <code>S_pq(t)=k(t,[T_p,T_q])</code> at discrete
    *  time <code>t</code> (continuous time <code>T_t</code>). 
    *  Needs all Libors <code>L_j, j&gt=p</code> up to time <code>t</code>.</p>
    *
    * <p>Straightforward implementation following the definition.
    * Extremely inefficient. Used only to test the correctness of the 
    * streamlined implementation.</p>
    *
    * @param p swap begins at <code>T_p</code>.
    * @param q swap ends at <code>T_q</code> (settled).
    * @param t current (discrete) time.
    */ 
     public double swRate(int p, int q, int t)
     { 
        double num,        // numerator in the swap rate definition
               denom;      // denominator
        
        num=B(p,t)-B(q,t);
        denom=0;
        for(int j=p;j<q;j++)denom+=delta[j]*B(j+1,t);
       
        return num/denom;
     } //end Swap_Rate 
     
     
   /** <p>The forward swap rate <code>S_pq(t)=k(t,[T_p,T_q])</code> at discrete 
    *  time <code>t</code> (continuous time <code>T_t</code>). 
    *  Needs all Libors <code>L_j, j&gt=p</code> up to time <code>t</code>.</p>
    *
    * @param p swap begins at <code>T_p</code>.
    * @param q swap ends at <code>T_q</code> (settled).
    * @param t current (discrete) time.
    */ 
     public double swapRate(int p, int q, int t)
     { 
        double f=1+X[q-1][t], S=delta[q-1];
        for(int k=q-2;k>=p;k--){ S+=delta[k]*f; f*=(1+X[k][t]); }
 
        return (f-1)/S;
     } //end Swap_Rate


   /** <p>The forward swap rate <code>S_pq(t)=k(t,[T_p,T_q])</code> at time 
    *  <code>t=0</code>.</p>
    *
    * @param p swap begins at <code>T_p</code>.
    * @param q swap ends at <code>T_q</code> (settled).
    */ 
     public double swapRate(int p, int q)
     { 
        double f=1+x[q-1], S=delta[q-1];
        for(int k=q-2;k>=p;k--){ S+=delta[k]*f; f*=(1+x[k]); }
 
        return (f-1)/S;
     } //end Swap_Rate

     
/*******************************************************************************
 *
 *                         ANNUITY NUMERAIRE
 *
 ******************************************************************************/                    
     
   /** <p>The annuity <code>B_pq(t)=sum_{k=p}^{q-1}delta_kB_{k+1}(t)</code>.
    *  Needs all Libors <code>L_j, j&gt=p+1</code> up to time <code>t</code>.
    *  </p>
    *
    * <p>Straightforward implementation following the definition.
    * Extremely inefficient. Used only to test the correctness of the 
    * streamlined implementation.</p>
    *
    * @param p annuity begins at <code>T_p</code>.
    * @param q annuity ends at <code>T_q</code> (settled).
    * @param t current (discrete) time.
    */ 
     public double b_pq(int p, int q, int t)
     { 
         double S=0;
         for(int k=p;k<q;k++)S+=delta[k]*B(k+1,t);
         return S;
     } //end B_pq
     
     
   /** <p>The annuity <code>B_pq(t)=sum_{k=p}^{q-1}delta_kB_{k+1}(t)</code>.
    *  Needs all Libors <code>L_j, j&gt=p+1</code> up to time <code>t</code>.
    *  </p>
    *
    * @param p annuity begins at <code>T_p</code>.
    * @param q annuity ends at <code>T_q</code> (settled).
    * @param t current (discrete) time.
    */ 
     public double B_pq(int p, int q, int t)
     { 
         double S=0, F=B(q,t);
         for(int k=q-1;k>=p;k--){ S+=delta[k]*F; F*=(1+X[k][t]); }

         return S;
     } //end B_pq

     
   /** <p>The annuity <code>B_pq(t)=sum_{k=p}^{q-1}delta_kB_{k+1}(t)</code>
    *  at time <code>t=0</code>.
    *  </p>
    *
    * @param p annuity begins at <code>T_p</code>.
    * @param q annuity ends at <code>T_q</code> (settled).
    */ 
     public double B_pq(int p, int q)
     { 
         double S=0, F=B0(q);
         for(int k=q-1;k>=p;k--){ S+=delta[k]*F; F*=(1+x[k]); }
         return S;
     } //end B_pq

        
/*******************************************************************************
 *
 *                         toString()
 *
 *******************************************************************************/        
    
    /** A message what type of factor loading it is, all the parameter values.
     */
    public String toString()
    {
        String stringRep=factorLoading.toString();
        stringRep+="\n\nInitial Libors: "+new ColtVector(l).toString();
        stringRep+="\n\nLibor volatilities:\n";
        
        for(int i=1;i<n;i++){
            
            double vol_i=vol(i);
            vol_i=FinMath.round(vol_i,3);
            stringRep+=vol_i+", ";
            if((i%20==0)||(i==n-1))stringRep+=vol_i+"\n";
        }
        
        return stringRep;
     }
             
/*******************************************************************************
 *
 *                           TEST PROGRAM
 *
 ******************************************************************************/            
    
    /** <p>Test program. Allocates sample Libor process in dimension 
     *  <code>n=20</code>
     *  then computes a histogram of <code>L_12(T_12)</code>.</p>
     *
     *  <p>Sample size: 50,000 paths.</p>
     */   
    public static void main(String[] args)
    {
        double before, after, time;
        before=System.currentTimeMillis();
        
        int n=20,            // dimension of Libor process
            Npath=50000,     // number of paths
            nBins=500,       // number of histogram bins
            i=12;            // Libor index
        
        LMM_Parameters lmmParams=new LMM_Parameters(n,LMM_Parameters.CS);
        LiborProcess LP=new LiborProcess(lmmParams);
        
        RandomVariable L_i=LP.L(i);
        String title="Forward Libor",
               axisLabel="L_"+i+"(T_"+i+")";
        L_i.displayHistogram(Npath,nBins,true,title,axisLabel);

        after=System.currentTimeMillis();
        time=after-before;
        System.out.println
        (Npath+" Paths of L_j(t), j="+i+",...,"+(n-1)+" up to time t=T_"+i+ 
         "\ntime: "+time);
        

    } // end main
    


} // end LiborProcess


/* WARNING: when the Libor random variable are programmed with path segments
 * the program crashes or does nor terminate. Investigate this.
 */