A few months ago I posted this pgplot powered fortran animation program on calculating the approximate value of PI by Monte Carlo method using random numbers.

And the other found myself reading about random or stochastic phenomenons. After reading about stochastic processes, I was invariably reminded of the stock market. The most stochastic place one can imagine.

After a bit of googling and a bit of tinkering, created the following program.

Given the annual drift, volatility, initial price and number of trials, the program simulates the movement of a stock price.

program MonteCarloStocks
double precision a,b,c
character(len=100) buffer
integer i,n,noarg,IARGC
noarg=IARGC()
if(noarg == 4) then
CALL GETARG (1, buffer)
read(buffer,'(f6.2)') a
CALL GETARG (2, buffer)
read(buffer,'(f6.2)') b
CALL GETARG (3, buffer)
read(buffer,'(f6.2)') c
CALL GETARG (4, buffer)
read(buffer,'(i7)') i
print *, a,b,c,i
else
write(*,'(a)',advance='no') "Please enter annual drift :"
read(*,'(f8.2)')a
write(*,'(a)',advance='no') "Please enter annual Volatility :"
read(*,'(f8.2)')b
write(*,'(a)',advance='no') "Please enter Initial Price :"
read(*,'(f8.2)')c
write(*,'(a)',advance='no') "Please enter no of iteration :"
read(*,'(i7)')i
end if
call monte(a,b,c,i)
call system("stock.xls")
end program
subroutine monte(driftyearly,volatility,initialprice,not)
double precision driftyearly,volatility,initialprice
double precision driftdaily,voldaily,p,pri(252),PPND16
integer,parameter :: nod=252
integer not
driftyearly = driftyearly/100
volatility = volatility/100
driftdaily=driftyearly/nod
voldaily=volatility/(nod**0.5)
open(200,file="fort.200")
call random_seed
if(not <= 0) not=1
if(not >=25000) not =25000
do j=1,not
pri(1)=initialprice
do i=2,nod
call random_number(p)
p=PPND16 (p, IFAULT)
pri(i)=pri(i-1)*EXP(driftdaily+voldaily*p)
end do
write(200,'(252f8.2)')pri
end do
close(200)
end subroutine
DOUBLE PRECISION FUNCTION PPND16 (P, IFAULT)
!
! ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3
!
! Produces the normal deviate Z corresponding to a given lower
! tail area of P; Z is accurate to about 1 part in 10**16.
!
! The hash sums below are the sums of the mantissas of the
! coefficients. They are included for use in checking
! transcription.
!
DOUBLE PRECISION ZERO, ONE, HALF, SPLIT1, SPLIT2, CONST1, &
CONST2, A0, A1, A2, A3, A4, A5, A6, A7, B1, B2, B3, &
B4, B5, B6, B7, &
C0, C1, C2, C3, C4, C5, C6, C7, D1, D2, D3, D4, D5, &
D6, D7, E0, E1, E2, E3, E4, E5, E6, E7, F1, F2, F3, &
F4, F5, F6, F7, P, Q, R
PARAMETER (ZERO = 0.D0, ONE = 1.D0, HALF = 0.5D0, &
SPLIT1 = 0.425D0, SPLIT2 = 5.D0, &
CONST1 = 0.180625D0, CONST2 = 1.6D0)
!
! Coefficients for P close to 0.5
!
PARAMETER (A0 = 3.3871328727963666080D0, &
A1 = 1.3314166789178437745D+2,&
A2=1.9715909503065514427D+3,&
A3=1.3731693765509461125D+4,&
A4=4.5921953931549871457D+4,&
A5=6.7265770927008700853D+4,&
A6=3.3430575583588128105D+4,&
A7=2.5090809287301226727D+3,&
B1=4.2313330701600911252D+1,&
B2=6.8718700749205790830D+2,&
B3=5.3941960214247511077D+3,&
B4=2.1213794301586595867D+4,&
B5=3.9307895800092710610D+4,&
B6=2.8729085735721942674D+4,&
B7=5.2264952788528545610D+3)
! HASH SUM AB 55.88319 28806 14901 4439
!
! Coefficients for P not close to 0, 0.5 or 1.
!
PARAMETER (C0 = 1.42343711074968357734D0,&
C1=4.63033784615654529590D0,&
C2=5.76949722146069140550D0,&
C3=3.64784832476320460504D0,&
C4=1.27045825245236838258D0,&
C5=2.41780725177450611770D-1,&
C6=2.27238449892691845833D-2,&
C7=7.74545014278341407640D-4,&
D1=2.05319162663775882187D0,&
D2=1.67638483018380384940D0,&
D3=6.89767334985100004550D-1,&
D4=1.48103976427480074590D-1,&
D5=1.51986665636164571966D-2,&
D6=5.47593808499534494600D-4,&
D7=1.05075007164441684324D-9)
! HASHSUMCD49.33206503301610289036
!
! CoefficientsforPnear0or1.
!
PARAMETER(E0=6.65790464350110377720D0,&
E1=5.46378491116411436990D0,&
E2=1.78482653991729133580D0,&
E3=2.96560571828504891230D-1,&
E4=2.65321895265761230930D-2,&
E5=1.24266094738807843860D-3,&
E6=2.71155556874348757815D-5,&
E7=2.01033439929228813265D-7,&
F1=5.99832206555887937690D-1,&
F2=1.36929880922735805310D-1,&
F3=1.48753612908506148525D-2,&
F4=7.86869131145613259100D-4,&
F5=1.84631831751005468180D-5,&
F6=1.42151175831644588870D-7,&
F7=2.04426310338993978564D-15)
! HASH SUM EF 47.52583 31754 92896 71629
!
IFAULT = 0
Q = P - HALF
IF (ABS(Q) .LE. SPLIT1) THEN
R = CONST1 - Q * Q
PPND16 = Q * (((((((A7 * R + A6) * R + A5) * R + A4) * R + A3) &
* R + A2) * R + A1) * R + A0) / &
(((((((B7 * R + B6) * R + B5) * R + B4) * R + B3) &
* R + B2) * R + B1) * R + ONE)
RETURN
ELSE
IF (Q .LT. ZERO) THEN
R = P
ELSE
R = ONE - P
END IF
IF (R .LE. ZERO) THEN
IFAULT = 1
PPND16 = ZERO
RETURN
END IF
R = SQRT(-LOG(R))
IF (R .LE. SPLIT2) THEN
R = R - CONST2
PPND16 = (((((((C7 * R + C6) * R + C5) * R + C4) * R + C3) &
* R + C2) * R + C1) * R + C0) / &
(((((((D7 * R + D6) * R + D5) * R + D4) * R + D3) &
* R + D2) * R + D1) * R + ONE)
ELSE
R = R - SPLIT2
PPND16 = (((((((E7 * R + E6) * R + E5) * R + E4) * R + E3) &
* R + E2) * R + E1) * R + E0) / &
(((((((F7 * R + F6) * R + F5) * R + F4) * R + F3) &
* R + F2) * R + F1) * R + ONE)
END IF
IF (Q .LT. ZERO) PPND16 = - PPND16
RETURN
END IF
END

Cool I would say.

An example of Monte Carlo method in action.

### Like this:

Like Loading...