% 1.731 event probabilities from Monte Carlo simulation
function runoff_monte
n=100
nrep=1000
threshold=1.5
% initial flow = 1
q(1)=0.5;
event_count=0;
% start Monte Carlo loop
for i=1:nrep
    % generate random rainfall for 100 days
    % prob(rain) on any given day = 0.2
    % generate precip values from -4 to 1
    precip=unifrnd(-20,5,1,n);
    % set all negative precip values to 0
    precip=precip.*(precip > 0);
    % start time loop
    for t=1:n
        q(t+1)=0.8*q(t)+0.2*precip(t);
    end
    qbar(i)=mean(q);
    % test to see if threshold exceeded 3 times during 100 day
    % simulation period ... if so, increment event counter
    if(sum(q > threshold)>3) event_count=event_count+1; end
end
qmean=mean(qbar)
prob_event=event_count/nrep
return