qmcintegration.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]

*/

/*
 * QmcIntegration.java
 *
 * Created on May 24, 2002, 3:22 PM
 */

package Examples.QuasiMonteCarlo;

import QuasiRandom.*;
import com.skylit.io.EasyReader;
import IO.*;
import Statistics.*;

/** <p>Tests the efficiency of the low discrepancy sequnces by computing
 *  the QMC integral of various {@link QuasiRandom.Intgrnd}s 
 *  over the unit cube Q=(0,1)^dim.
 *  The relative error is reported as a signed prcentage of the analytic value
 *  of the integral. <i>The analytic integral must be nonzero</i>.</p>
 *
 * <p>The QMC integrals are evaluated for N=2^n-1 low discrepancy points.
 * At these values of N the Gray code counter "is in sync" with the ordinary
 * sequence 1,2,3,...,N again:<br>
 * {1,2,...,N}={Gray(1),Gray(2),...,Gray(N)}.</p>
 *
 * <p>The results are written to a file "QmcIntDim_dim" in the current 
 * directory.</p>
 *
 * @author  Michael J. Meyer
 */
public class QmcIntegration {
 
   static int dim;

   
   public static void main(String[] args)
   {       
        EasyReader io=new EasyReader();

        System.out.println
        ("TESTING THE LOW DISCREPANCY SEQUENCES.\n\n"+
         "Enter dimension: ");

        dim=io.readInt();
        
        // file writer to write in columns of fixed width w
        final int w=14;                     // width of output columns
        String fileName="QmcIntDim"+dim;    // output file
        FixedFieldWidthFileWriter fout=null;
        
        try{ fout=new FixedFieldWidthFileWriter(fileName,w); }
        catch(java.io.IOException e)
        { System.out.println("Cannot allocate file "+fileName+", exiting...");
          System.exit(0); }
            
         String message=
         "QUASI MONTE CARLO COMPUTATION OF TEST INTEGRALS:\n\n"+
         "We compute Int_Q F(x)dx, where Q=(0,1)^d (unit cube), d=dimension\n"+
         "and F(x)=h(x_1)h(x_2)...h(x_d), where h(u)=g(m*u-[m*u]).\n\n"+
         "Here [t] is the greatest integer <= t as usual. The function h(u)\n"+
         "is periodic with m periods on (0,1) repeating the function g(u)\n"+
         "m times on (0,1). If g(u) has one peak on (0,1), h(u) has m peaks\n"+
         "and F(x) has m^d (very narrow) peaks on Q.\n\n"+
         "DIMENSION="+dim;
             
   try{ 
            
        fout.write(message);
        
        // the integrands
        java.util.Vector integrands=new java.util.Vector();
        
        integrands.add(new Intgrnd_1(dim,1));
        integrands.add(new Intgrnd_1(dim,2));
        integrands.add(new Intgrnd_1(dim,3));
        
        integrands.add(new Intgrnd_2(dim,1));
        integrands.add(new Intgrnd_2(dim,2));
        integrands.add(new Intgrnd_2(dim,3));
        
        integrands.add(new Intgrnd_3(dim,1));
        integrands.add(new Intgrnd_3(dim,2));
        integrands.add(new Intgrnd_3(dim,3));
        
        integrands.add(new Intgrnd_4(dim,1,2));
        integrands.add(new Intgrnd_4(dim,2,3));
        integrands.add(new Intgrnd_4(dim,3,4));
        
        integrands.add(new Intgrnd_5(dim,1));
        integrands.add(new Intgrnd_5(dim,2));
        integrands.add(new Intgrnd_5(dim,3));

        integrands.add(new Intgrnd_6(dim,1));
        integrands.add(new Intgrnd_6(dim,2));
        integrands.add(new Intgrnd_6(dim,3));
       
        
        for(int i=0;i<integrands.size();i++)
        {
            CubeFunction F=(CubeFunction)integrands.elementAt(i);
            fout.write("\n\n\n"+F.getName()+", QMC relative error (%):\n\n");
            computeIntegrals(F,fout); 
        }
        
        fout.close();
        System.out.println
        ("Finished, results in file "+fileName+" in current directory.");
 
   }// end try
   catch(java.io.IOException e)
   { System.out.println("Cannot write to file "+fileName+", exiting...");
     System.exit(0); }

     } // end main


     // QMC computation of Integral_Q F(x)dx using the Sobol sequence S,
     // the Niederreiter-Xing sequence NXS, the Halton sequence and a 
     // MersenneTwister uniform sequence U.
     public static void 
     computeIntegrals(CubeFunction F, FixedFieldWidthFileWriter fout)
     throws java.io.IOException
     {
        int N=32767;         // number of points N=2^15-1
        
        LowDiscrepancySequence SS=new Sobol(dim), 
                               NXS=new NX(dim), 
                               US=new Uniform(dim), 
                               HS=new Halton(dim);
         
        double ssum=0, nxssum=0, usum=0, hsum=0,
               sint, nxsint, uint, hint,             // QMC integrals
               serr, nxserr, uerr, herr;             // relative errors
        
        double aint=F.integral();    // analytic integral
        int twopower=512;            // the smallest power of two >=n
        

        for(int n=0;n<N;n++)       // n=2^k
        {
           ssum+=F.value(SS.nextPoint());
           nxssum+=F.value(NXS.nextPoint());
           hsum+=F.value(HS.nextPoint());
           usum+=F.value(US.nextPoint());
           
           if(n==(twopower-2))       // evaluate integrals and relative error
           {
              sint=FinMath.round(ssum/n,5);
              nxsint=FinMath.round(nxssum/n,5);
              hint=FinMath.round(hsum/n,5);
              uint=FinMath.round(usum/n,5);
              
              fout.writeField("N: "+(n+1));
              serr=FinMath.round(100*(sint-aint)/aint,2);
              fout.writeField("SOB: "+serr);
              nxserr=FinMath.round(100*(nxsint-aint)/aint,2);
              fout.writeField("NX: "+nxserr);
              herr=FinMath.round(100*(hint-aint)/aint,2);
              fout.writeField("HAl: "+herr);
              uerr=FinMath.round(100*(uint-aint)/aint,2);
              fout.write("UNF: "+uerr+"\n");
              
              twopower*=2;    // set this to the next power of two.
           }
         } // end for n
         
        
      } // end computeIntegrals
     

     // control output field width
     public static void write(String str, int fieldwidth)
     {
         System.out.print(str);
         int blanks=fieldwidth-str.length();
         for(int j=0;j<blanks;j++)System.out.print(" ");
     }
     


}// end QmcIntegration