deltaT = 0.25 ## of a year
sDeltaT = sqrt(deltaT)

g0 = 400 
r0 = 0.04

n_simulations = 10000
n_steps = 10/deltaT 

## Simulate paths of g and r and then S:
##
all_gs = matrix(data=rep(NA, n_simulations*(n_steps+1)), nrow=n_simulations, ncol=(n_steps+1))
all_rs = matrix(data=rep(NA, n_simulations*(n_steps+1)), nrow=n_simulations, ncol=(n_steps+1))
all_Ss = matrix(data=rep(NA, n_simulations*(n_steps+1)), nrow=n_simulations, ncol=(n_steps+1))
for( ii in 1:n_simulations ){

    ## Forward sweeps:
    ## 
    gs = rep(NA, n_steps+1); gs[1] = g0
    rs = rep(NA, n_steps+1); rs[1] = r0
    for( ti in 1:n_steps ){
        dg = rs[ti] * gs[ti] * deltaT + 0.25 * gs[ti] * sDeltaT * rnorm(1)
	dr =  0.005 * rs[ti] * deltaT + 0.01 * rs[ti] * sDeltaT * rnorm(1)
	gs[ti+1] = gs[ti] + dg
	rs[ti+1] = rs[ti] + dr
    }

    ## Backwards sweeps:
    ## 
    Ss = rep(0, n_steps+1)
    for( ti in seq(n_steps, 1, by=-1) ){
        c = max(gs[ti+1] - 200, 0) * 10000
	dS = rs[ti+1] * Ss[ti+1] * deltaT + c * deltaT
	Ss[ti] = Ss[ti+1] + dS
    }

    ## Save things:
    ## 
    all_gs[ii, ] = gs
    all_rs[ii, ] = rs
    all_Ss[ii, ] = Ss    
}

if( F ){
    ##
    ## This is fixed with the Randleman and Bartter model: 
    ##
    negative_rate_rows = apply(all_rs, 1, function(row) any(row < 0)) ## I think negative rates are not valid so I drop them
    all_rs = all_rs[!negative_rate_rows, ]
    all_gs = all_gs[!negative_rate_rows, ]
    all_Ss = all_Ss[!negative_rate_rows, ]
}

hist(all_Ss[, 1]/1e6)
print(sprintf('n_simulations= %d; estimated Continuco gold mine= %10.6f', n_simulations, mean(all_Ss[, 1]/1e6)))