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

*/




/*
 * DirichletProblem.java
 *
 * Created on January 22, 2002, 8:00 PM
 */
 
 
package Examples.Probability;

import Statistics.*;
import Processes.*;
import java.lang.Math.*;


 

/** <p> We solve the Dirichlet problem on the unit disk with boundary function 
 * h(x,y)=x*y (which is harmonic on the entire plane). Thus the solution 
 * is given by the same formula on the interior of the disc.</p>
 * 
 * <p>Fixing the point x=(1/4,1/4) we compute the solution f(x) as
 * f(x)=E(h(B_tau)), where B is a two dimensional Brownian motion 
 * started at the point x and tau the first time B hits the boundary 
 * of the disc.</p>
 *
 * <p>As a proxy for the point at which the Brownian motion hits the boundary
 * we take the first point a past the boundary projected radially back 
 * on the circle, ie.: a/||a||.</p>
 *
 * @author  Michael J. Meyer
 */
public class DirichletProblem{

    /** projects the point (u,v) radially on the unit circle */
    public static double[] boundaryProjection(double[] x)
    {
        double u=x[0], v=x[1], f=Math.sqrt(u*u+v*v);
        double[] result={u/f,v/f};
        return result;
    } 
     
     
    public static void main (String[] args)
    {
        
        int T=50000,      //time steps to horizon
            dim=2;        //dimension
        double dt=0.001;  //size of time step
        
        double[] x={0.25,0.25};  //starting point
        
        
        //allocate a two dimensional Brownian motion starting at x
        final VectorBrownianMotion B=new VectorBrownianMotion(dim,T,dt,x);
        
        
        //define the unit disc as a two dimensional region
        Region_nD disc=new Region_nD(){
        
             public boolean isMember(double[] z)
             {
                  double u=z[0], v=z[1];
                  return (u*u+v*v<1);       }
                  
        }; //end disc
        
        
        //allocate the hitting time for the boundary
        final StoppingTime tau=new FirstExitTime_nD(B,disc);
        
        
        // allocate the random variable h(B_tau) 
        RandomVariable hB_tau=new RandomVariable(){
        
            public double getValue(int t)
            {
                double[] z=boundaryProjection(B.sampledAt(tau).getValue(t));
                double u=z[0], v=z[1];
                return u*v;
             }
             
         }; // end hB_tau
        
        
        
        //compute f(x)=E(h(B_tau)) over 20000 paths 
        //this should be 0.25*0.25=0.0625
        double fx=hB_tau.expectation(40000);
        //cut this down to 5 decimals
        fx=1.0*Math.round(10000*fx)/10000;
        
        String message="The solution f(x) at x=(1/2,1/2) should be 0.0625\n"+
                       "and is computed as "+fx;
                       
        System.out.println(message);
        

    } // end main
    
} // end DirichletProblem


/** REMARK:
 *  The horizon T has to be chosen large enough and the time step dt small enough 
 *  for this to work well. The computation does take some time!
 */