Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members  

StochasticProcess.h

Go to the documentation of this file.
00001 /* WARANTY NOTICE AND COPYRIGHT
00002 This program is free software; you can redistribute it and/or
00003 modify it under the terms of the GNU General Public License
00004 as published by the Free Software Foundation; either version 2
00005 of the License, or (at your option) any later version.
00006 
00007 This program is distributed in the hope that it will be useful,
00008 but WITHOUT ANY WARRANTY; without even the implied warranty of
00009 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00010 GNU General Public License for more details.
00011 
00012 You should have received a copy of the GNU General Public License
00013 along with this program; if not, write to the Free Software
00014 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
00015 
00016 Copyright (C) Michael J. Meyer
00017 
00018 matmjm@mindspring.com
00019 spyqqqdia@yahoo.com
00020 
00021 */
00022 
00023 #ifndef martingale_stochasticprocess_h
00024 #define martingale_stochasticprocess_h
00025 
00026 // template class fully defined in header
00027 #include "Array.h"
00028 #include "RandomObject.h"
00029 #include "StochasticGenerator.h"
00030 
00036 // StoppingTime, Region, HittingTime, FirstExitTime, StochasticProcess, PathFunctional. 
00037 
00038 
00039 MTGL_BEGIN_NAMESPACE(Martingale)
00040 
00041 
00042 
00043 /********************************************************************************************
00044  *
00045  *                STOPPING TIMES, HITTING TIMES, FIRST EXIT TIMES
00046  * 
00047  ********************************************************************************************/
00048  
00049 
00051 class StoppingTime {
00052         
00053 protected:
00054         
00055         int T;  // number of time steps to horizon
00056         
00057 public:
00058         
00060         StoppingTime(int T_oo) : T(T_oo) {  }
00061 
00066     virtual bool stop(int t) = 0;
00067         
00068    
00071    int valueAlongCurrentPath()
00072    {
00073            int t=0; while(!stop(t)&&(t<T)) t++;
00074            return t;
00075    }
00076     
00077 }; //end StoppingTime
00078 
00079 
00080 
00083 template<typename RangeType>
00084 class Region {
00085         
00086 public:
00087 
00090   virtual bool isMember(const RangeType& x) = 0;
00091   
00092   
00093 }; //end Region
00094 
00095 
00096 template<typename RangeType, typename ScalarType>
00097 class StochasticProcess;     // defined below
00098 
00099 
00107 template<typename RangeType, typename ScalarType>
00108 class HittingTime : public StoppingTime {
00109 
00110    StochasticProcess<RangeType,ScalarType>* X;
00111    Region<RangeType>* D;
00112         
00113 public:
00114    
00118    HittingTime
00119    (StochasticProcess<RangeType,ScalarType>* Y, Region<RangeType>* G) : 
00120    StoppingTime(Y->getT()), X(Y), D(G) 
00121    {  }
00122 
00123    
00125    bool stop(int t){ return ((D->isMember(X->currentPath(t)))||(t==T)); }
00126           
00127 }; //end HittingTime
00128 
00129 
00130 
00131 
00139 template<typename RangeType, typename ScalarType>
00140 class FirstExitTime : public StoppingTime {
00141 
00142    StochasticProcess<RangeType,ScalarType>* X;
00143    Region<RangeType>* D;
00144         
00145 public:
00146    
00150    FirstExitTime
00151    (StochasticProcess<RangeType,ScalarType>* Y, Region<RangeType>* G) : 
00152    StoppingTime(Y->getT()), X(Y), D(G) 
00153    {  }
00154 
00155    
00157    bool stop(int t){ return (!(D->isMember(X->currentPath(t)))||(t==T)); }
00158        
00159        
00160 }; //end FirstExitTime
00161 
00162 
00163 
00164 /********************************************************************************************
00165  *
00166  *                       REGIONS IN EUCLIDEAN SPACE
00167  * 
00168  ********************************************************************************************/ 
00169 
00172 class EuclideanRegion : public Region< RealVector > {
00173         
00174 protected:
00175         
00176         int dim;  // dimension
00177         
00178 public:
00179         
00180         int getDimension() const { return dim; }
00181         
00183         EuclideanRegion(int d) : dim(d) { }
00184         
00186         virtual RealVector& boundaryProjection(const RealVector& x) = 0;
00187         
00193         RealVector& boundaryIntersection
00194         (const RealVector& u, const RealVector& v, int N=5) 
00195     {
00196                 RealVector x(u), y(v), z(u);
00197                 for(int i=0;i<N;i++){          
00198                         
00199                         z=x; z+=y; z*=0.5;
00200                         if(isMember(z)) x=z; else y=z;
00201                 }
00202                 return boundaryProjection(y);
00203         } // end boundaryIntersection
00204                         
00205 }; // end EuclideanRegion
00206 
00207 
00208 // BALL CENTERED AT THE ORIGIN
00209 
00212 class Ball : public EuclideanRegion {
00213         
00214         Real R;         // radius
00215         
00216 public:
00217         
00218         Real getRadius() const { return R; }
00219         
00223         Ball(int d, Real r=1.0) : EuclideanRegion(d), R(r) {  }
00224         
00225         bool isMember(const RealVector& x)
00226     {
00227                 Real sum=0;
00228                 for(int i=0;i<dim;i++) sum+=x[i]*x[i];
00229                 return(sum<R*R);
00230         }
00231         
00233         RealVector& boundaryProjection(const RealVector& x)
00234     {
00235                 RealVector& u=*(new RealVector(x));
00236                 u*=R/x.norm();
00237                 return u;
00238         }
00239         
00240 }; // end Ball
00241          
00242         
00243         
00244         
00245 
00246 
00247 /********************************************************************************************
00248  *
00249  *                            STOCHASTIC PROCESS
00250  * 
00251  ********************************************************************************************/
00252  
00253 
00286 template<typename RangeType=Real, typename ScalarType=Real>
00287 class StochasticProcess {
00288         
00289 protected:
00290         
00291         int dim;      // dimension of range
00292     int T;        // number of time steps to horizon
00293         
00294 public:
00295         
00297     int getDimension() const { return dim; }
00298     
00300     int getT() const { return T; }
00301         
00306     virtual RangeType& currentPath(int t) = 0;
00307  
00308         
00309 // CONSTRUCTOR
00310         
00314     StochasticProcess(int T_oo) : dim(1), T(T_oo) {  }
00315     
00320     StochasticProcess(int d, int T_oo) : dim(d), T(T_oo) {  }
00321         
00322         virtual ~StochasticProcess(){ }
00323 
00324                                              
00325 // TIME STEPS
00326         
00332     virtual void timeStep(int t) = 0;
00333         
00334    
00335 
00350     virtual void timeStep(int t, int s){  for(int u=t;u<s;u++)timeStep(u);  }
00351     
00352 
00353 // PATHS
00354 
00363     virtual void newPathSegment(int t, int s){  for(int u=t;u<s;u++)timeStep(u);  }
00364     
00365       
00366     
00374     virtual int newPathSegment(int t, StoppingTime* tau)
00375     {
00376         int s=t;
00377         while(!(tau->stop(s))){ timeStep(s); s++; }
00378         return s;
00379      }
00380      
00381      
00387     virtual void newPathSegment(int t){  newPathSegment(0,t);  }
00388         
00389     
00396     virtual int newPathSegment(StoppingTime* tau){  return newPathSegment(0,tau);  }
00397         
00398         
00409     virtual void newPathBranch(int t){ newPathSegment(t,T); }
00410         
00411 
00414     virtual void newPath(){ newPathBranch(0); }
00415     
00416 
00417 // SAMPLING
00418         
00426     RandomObject<RangeType,ScalarType>* operator()(int t, StoppingTime* tau)
00427     {
00428         return new SampledAt(this,t,tau);
00429         } 
00430     
00436     RandomObject<RangeType,ScalarType>* operator()(StoppingTime* tau)
00437     {
00438         return new SampledAt(this,t,tau);
00439         } 
00440         
00441 private:
00442         
00447         class SampledAt : public RandomObject<RangeType,ScalarType> {
00448                 
00449                 StochasticProcess X;
00450                 StoppingTime tau;
00451                 int t;
00452                 
00453          public:
00454                 
00455                 SampledAt(StochasticProcess* Y, int s, StoppingTime* psi) : 
00456                 RandomObject<RangeType,ScalarType>(Y->getDimension()), X(Y), tau(psi), t(s) 
00457                 {  }
00458                 RangeType nextValue(){  return X->currentPath(X->newPathSegment(t,tau)); }
00459                 
00460         }; // end SampledAt
00461     
00462     
00463 }; //end StochasticProcess
00464 
00465 
00466 
00467 /********************************************************************************************
00468  *
00469  *                    VECTOR VALUED STOCHASTIC PROCESS
00470  * 
00471  ********************************************************************************************/
00472 
00473 typedef StochasticProcess< RealVector, Real > VectorProcess;
00474 
00475 
00496 class BrownianVectorProcess : public VectorProcess {
00497 
00498 protected:
00499         
00500         // rows addressed as Real*
00501         Array1D< RealVector* > path;      // *(path[t]) is the state of the process at time t.
00502         RealMatrix Z;                     // the row Z[t][] is the vector driving the time step t->t+1. 
00503         StochasticGenerator* SG;          // the generator computing the Z-matrix
00504         
00505 public:
00506         
00508         RealVector& currentPath(int t) { return *(path[t]); }
00509 
00511         void setStart(const RealVector& x){ currentPath(0)=x; }
00512  
00513 // CONSTRUCTOR
00514         
00518         BrownianVectorProcess(int dim, int T) : VectorProcess(dim,T),
00519         path(T+1), Z(T,dim,0,0), 
00520         SG(new MonteCarloVectorDriver(dim))
00521     {   
00522                 for(int t=0;t<=T;t++) path[t]=new RealVector(dim);
00523         }
00524         
00525         ~BrownianVectorProcess()
00526         {  
00527                 for(int t=0;t<=T;t++) delete path[t];
00528                 delete SG; 
00529         }
00530 
00531 
00532 // PATH GENERATION
00533                 
00534 
00535     virtual void newPathSegment(int t, int s)
00536         {  
00537                 SG->newWienerIncrements(t,s,Z);
00538                 for(int u=t;u<s;u++)timeStep(u);  
00539         }
00540         
00546     virtual void newPathSegment(int t){  newPathSegment(0,t);  }
00547     
00548       
00549     virtual int newPathSegment(int t, StoppingTime* tau)
00550     {
00551         SG->newWienerIncrements(t,T,Z);
00552                 int s=t;
00553         while(!(tau->stop(s))){ timeStep(s); s++; }
00554         return s;
00555      }
00556          
00563     virtual int newPathSegment(StoppingTime* tau){  return newPathSegment(0,tau);  }
00564         
00567         void printCurrentPath(int t)
00568     {
00569                 cout << "\n\nBrownian vector process, current path:\n";
00570                 for(int s=0;s<=t;s++)
00571                 cout << "\nt = " << s << ", X(t):" << currentPath(s);
00572         }
00573         
00574          
00575 // switching between Monte Carlo and Quasi Monte Carlo dynamics
00576          
00579     void switchToQMC() 
00580         {  
00581                 if(SG) delete SG;
00582                 SG = new SobolVectorDriver(dim,T);
00583         }
00584         
00585         
00588     void switchToMC() 
00589         { 
00590                 if(SG) delete SG;
00591                 SG = new MonteCarloVectorDriver(dim);
00592         }
00593          
00594          
00595  }; // BrownianVectorProcess
00596      
00597     
00598    
00599 
00600 
00601 /********************************************************************************************
00602  *
00603  *                    REAL STOCHASTIC PROCESS
00604  * 
00605  ********************************************************************************************/
00606 
00607                 
00608 typedef StochasticProcess<Real,Real> ScalarProcess;
00609 
00610 
00616 class BrownianScalarProcess : public ScalarProcess {
00617         
00618 protected:
00619         
00620         RealVector path;          // path[t] is he state of the process at time t
00621         RealVector  Z;            // the row Z[t] is the standard normal increment driving 
00622                                   // the time step t->t+1.
00623         StochasticGenerator* SG;  // the generator computing the Z-vector
00624         
00625 public:
00626         
00628         RealVector& getPath() { return path; }
00629                         
00631         Real& currentPath(int t) { return path[t]; }
00632         
00633 // Constructor
00634         
00637         BrownianScalarProcess(int T) : ScalarProcess(T),
00638         path(T+1),
00639         Z(T), 
00640         SG(new MonteCarloScalarDriver())
00641     {   }
00642         
00643         ~BrownianScalarProcess(){ delete SG; }
00644 
00645 // the new path generation methods
00646         
00647 
00648     virtual void newPathSegment(int t, int s)
00649         {  
00650                 SG->newWienerIncrements(t,s,Z);
00651                 for(int u=t;u<s;u++)timeStep(u);  
00652         }
00653     
00654       
00655     virtual int newPathSegment(int t, StoppingTime* tau)
00656     {
00657         SG->newWienerIncrements(t,T,Z);
00658                 int s=t;
00659         while(!(tau->stop(s))){ timeStep(s); s++; }
00660         return s;
00661      }
00662          
00663 // switching between Monte Carlo and Quasi Monte Carlo dynamics
00664          
00667     void switchToQMC() 
00668         {  
00669                 if(SG) delete SG;
00670                 SG = new SobolScalarDriver(T);
00671         }
00672         
00673         
00676     void switchToMC() 
00677         { 
00678                 if(SG) delete SG;
00679                 SG = new MonteCarloScalarDriver();
00680         }
00681          
00682  }; // BrownianScalarProcess
00683      
00684     
00685    
00686                 
00687 
00688 
00689 /********************************************************************************************
00690  *
00691  *                              PATH FUNCTIONALS
00692  * 
00693  ********************************************************************************************/
00694  
00733 template<typename ProcessRangeType=Real,         
00734          typename RangeType=Real, 
00735                  typename ScalarType=Real>
00736 class PathFunctional {
00737          
00738 protected:
00739          
00740          // the process of which it is a functional
00741          StochasticProcess<ProcessRangeType,ScalarType>* X;        
00742          int dim;      // dimension of the RangeType of the functional.
00743          
00744 public: 
00745         
00747          int getDimension(){ return dim; }
00748          
00750          StochasticProcess<ProcessRangeType,ScalarType>* getProcess(){ return X; }
00751          
00755          PathFunctional
00756          (StochasticProcess<ProcessRangeType,ScalarType>* Y, int d=1) : 
00757          X(Y), dim(d) {  }
00758          
00761          virtual RangeType valueAlongCurrentPath() = 0;
00762          
00763          
00772          RandomObject<RangeType,ScalarType>* 
00773          conditionedAtTime(int t){ return new H_t(this,t); }
00774          
00782          RangeType conditionalExpectation
00783          (int t, int nPath, bool progressReport=false, string message=" ")
00784      {
00785                  if(progressReport)
00786                  return conditionedAtTime(t)->expectation(nPath,message);
00787                  // otherwise
00788                  return conditionedAtTime(t)->expectation(nPath);
00789          }
00790 
00791          
00792          
00799          RangeType* conditionalMeanAndVariance
00800          (int t, int nPath, bool progressReport=false, string message=" ")
00801      {
00802                  if(progressReport)
00803                  return conditionedAtTime(t)->meanAndVariance(nPath,message);
00804                  // otherwise
00805                  return conditionedAtTime(t)->meanAndVariance(nPath);
00806          }
00807          
00808          
00809          
00810 private: 
00811          
00816          class H_t : public RandomObject<RangeType> {
00817                  
00818                  PathFunctional<ProcessRangeType,RangeType>* H;
00819                  int t;                      // time of conditioning
00820                  
00821          public:
00822                  
00826                  H_t(PathFunctional<ProcessRangeType,RangeType>* G, int s=0) : 
00827                  RandomObject<RangeType>(), H(G), t(s)
00828              {   }
00829                  
00833                  RangeType nextValue()
00834              {
00835                          StochasticProcess<ProcessRangeType,ScalarType>* X=H->getProcess();
00836                          X->newPathBranch(t);
00837                          return H->valueAlongCurrentPath();
00838                  }
00839          }; // end H_t
00840          
00841  }; // end PathFunctional
00842  
00843  
00848 class SumFunctional : public PathFunctional< RealVector > {
00849                          
00850         VectorProcess* X; int T;                
00851 public:
00852         SumFunctional(VectorProcess* Y, int T_oo) : 
00853         PathFunctional< RealVector >(Y), X(Y), T(T_oo) {  }
00854                          
00855         Real valueAlongCurrentPath()
00856     {
00857             RealVector& x=X->currentPath(T);
00858                 return x[0]+x[1];
00859         }
00860 }; // end SumFunctional
00861 
00862 MTGL_END_NAMESPACE(Martingale)
00863 
00864 #endif
00865 

Generated on Mon Sep 22 02:16:32 2003 for Libor-Library by doxygen1.3-rc3