1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
| c-----------------------------------------------------------------------
SUBROUTINE STEP ( n, mudt, chol, sqdt, dwork, step )
c-----------------------------------------------------------------------
c
c INPUT
c n : number of elements integer
c mudt[n] : mu(i)*deltat, (n) double
c chol[n,n] : choleski factor of covariance (n*n) double
c sqdt : sqrt of dt (discret. step in years) double
c
c WORKSPACE
c dwork : ( 2*n ) double
c
c OUTPUT
c step[n] : generated step (n) double
c
c CALL
c PMLV : M*V = vector (M(n*n), V vector(n),
c gives M*V vector(n))
c SVVX : V1 + V2*X, V1 and V2 vectors(n), X scalar, give vector(n)
c
c-----------------------------------------------------------------------
c
IMPLICIT NONE
c
c arguments i/o
INTEGER n
DOUBLE PRECISION sqdt
DOUBLE PRECISION mudt(*), chol(*), step(*)
c
c workspaces
DOUBLE PRECISION dwork(*)
c
c local variables
INTEGER i, pdb, pdc
DOUBLE PRECISION mean, std
PARAMETER ( mean = 0.0d0, std = 1.0d0 )
c
c-----------------------------------------------------------------------
c
c initialization
c
c pointers for double precision work space : dwork
c -------------------------------------------------
pdb = 1
c pdb : pointer for noise vector(nelem), so (n) more,
pdc = pdb + ( n )
c pdc : pointer for work vector(nelem), so (n) more,
c
c Total size of dwork array = (n) + (n) = ( 2*n )
c
c-----------------------------------------------------------------------
c
c generate n-Gaussian number N(m,std)
CALL IVXRN ( n, mean, std, dwork(pdb) )
c
c sigma*dW, vector(n)
CALL PMLV ( n, chol, dwork(pdb), dwork(pdc) )
c
c [mu*dt] + [sigma*dW*sqrt(t)], vector(n)
CALL SVVX ( n, mudt, dwork(pdc), sqdt, step )
c
RETURN
END |
Partager