Simulating stock market drawdowns

With the U.S. stock market down about 18% year to date and having fallen 7 weeks in a row, financial publications are asking
(1) Is it a good time to invest in stocks?
(2) Is this the low point of stock prices?

It may be surprising that even if the answer to (1) is “yes”, the answer to (2) is “probably not”. One can simulate a random walk with plausible parameters for the mean and standard deviation of daily stock market returns, assuming normality, and compute the probability that on any day in the next year (about 252 trading days) stocks will close below the current level. For example, if one assumes 10% annualized returns with 16% standard deviation, termed “volatility”, (fairly optimistic parameters), there is only a 7% chance that the stock market will not close lower on some day over the next year.

         #paths         #steps         thresh      ann_drift ann_volatility
        1000000            252           0.00          10.08          15.87

    first_lower    probability          cumul
              0       0.068836       0.000000
              1       0.483941       0.483941
              2       0.121955       0.605896
              3       0.060490       0.666386
              4       0.037821       0.704207
              5       0.026629       0.730836
              6       0.019807       0.750643
              7       0.015581       0.766224
              8       0.012829       0.779053
              9       0.010593       0.789646
             10       0.009099       0.798745
             21       0.002974       0.851031
             42       0.001052       0.885094
             63       0.000556       0.899898
             84       0.000323       0.908613
            105       0.000238       0.914455
            126       0.000190       0.918738
            147       0.000109       0.922032
            168       0.000115       0.924550
            189       0.000083       0.926623
            210       0.000082       0.928356
            231       0.000071       0.929864
            252       0.000054       0.931164

Even with the benign return and volatility assumptions used, there is a 31% chance that the stock market will be 10% lower than the starting level at some point.

         #paths         #steps         thresh      ann_drift ann_volatility
        1000000            252         -10.00          10.08          15.87

    first_lower    probability          cumul
              0       0.694206       0.000000
              1       0.000000       0.000000
              2       0.000000       0.000000
              3       0.000000       0.000000
              4       0.000001       0.000001
              5       0.000004       0.000005
              6       0.000014       0.000019
              7       0.000036       0.000055
              8       0.000107       0.000162
              9       0.000197       0.000359
             10       0.000255       0.000614
             21       0.001942       0.014034
             42       0.002614       0.066781
             63       0.002103       0.116601
             84       0.001757       0.156750
            105       0.001398       0.189179
            126       0.001096       0.215487
            147       0.000885       0.237243
            168       0.000795       0.255386
            189       0.000683       0.270854
            210       0.000542       0.284038
            231       0.000497       0.295660
            252       0.000439       0.305794

Here is the code, which runs in a few seconds for 10^6 paths.

module bottom_mod
implicit none
integer, parameter :: dp = kind(1.0d0)
private
public :: dp, first_lower
contains
function first_lower(max_steps, thresh, drift, sd) result(ifirst)
! return the first time that the cumulative sum of normal variates
! with given mean (drift) and standard deviation (sd) is below thresh
integer      , intent(in) :: max_steps
real(kind=dp), intent(in) :: thresh
real(kind=dp), intent(in) :: drift
real(kind=dp), intent(in) :: sd
integer                   :: ifirst
integer                   :: i
real(kind=dp)             :: x
x = 0.0_dp
do i=1,max_steps
   x = x + drift + sd*random_normal()
   if (x < thresh) then
      ifirst = i
      return
   end if
end do
ifirst = 0 ! threshold never breached
end function first_lower
!
function random_normal() result(fn_val)
!   Generate a random normal deviate using the polar method. (Alan Miller)
!   Reference: Marsaglia,G. & Bray,T.A. 'A convenient method for generating
!              normal variables',Siam Rev.,vol.6,260-264,1964.
real(kind=dp)  :: fn_val
! Local variables
real(kind=dp)           :: u,sum
real(kind=dp),save      :: v,sln
logical,save   :: second = .false.
real(kind=dp),parameter :: one = 1.0d0,vsmall = tiny(one)
if (second) then
! If second,use the second random number generated on last call
  second = .false.
  fn_val = v*sln
else
! First call; generate a pair of random normals
  second = .true.
  do
    call random_number(u)
    call random_number(v)
    u = scale(u,1) - one
    v = scale(v,1) - one
    sum = u*u + v*v + vsmall         ! vsmall added to prevent log(zero) / zero
    if(sum < one) exit
  end do
  sln = sqrt(- scale(log(sum),1) / sum)
  fn_val = u*sln
end if
end function random_normal
end module bottom_mod
program xbottom
use bottom_mod, only: dp, first_lower
implicit none
integer, parameter       :: npaths = 10**6, max_steps = 252
integer                  :: i, ipath
integer, allocatable     :: ifirst(:)
real(kind=dp), parameter :: trading_days_year = 252.0_dp, xscale = 100.0_dp, &
   drift = 0.0004_dp, sd = 0.01_dp, thresh = 0.0_dp
real(kind=dp) :: prob, cumul
allocate (ifirst(npaths))
print "(*(a15))","#paths","#steps","thresh","ann_drift","ann_volatility"
print "(2i15,*(f15.2))",npaths,max_steps, &
      xscale*[thresh,trading_days_year*drift,sqrt(trading_days_year)*sd]
do ipath=1,npaths
   ! for each path compute the first day that stocks closed below threshold
   ifirst(ipath) = first_lower(max_steps, thresh, drift, sd)
end do
print "(/,*(a15))","first_lower","probability","cumul"
do i=0,max_steps
   prob = count(ifirst==i)/real(npaths,kind=dp)
   if (i > 0) cumul = cumul + prob ! cumul prob. that threshold breached by day i
   if (i < 11 .or. mod(i,21) == 0) then
      print "(i15,2f15.6)",i,prob,cumul
   else if (i == 0) then 
      print "(i15,2f15.6)",i,prob ! probability that threshold never breached
   end if
end do
end program xbottom
4 Likes