% Matlab(R) code for 
% "day-to-day model" 
% as presented in 
% "A statistically predictive model for future monsoon failure in India"
% by J. Schewe and A. Levermann
%
% Author:	Jacob Schewe 
% Address:  Potsdam Institute for Climate Impact Research, Potsdam, Germany
% email:	jacob.schewe@pik-potsdam.de
% Date:		03 Jun 2012
%
% Please refer to the manuscript for a discussion of the model. 
%
clear all

% --- Setup ---
% --- Parameters
nyears = 6030; % Number of years (seasons) to run
lseason = 135; % Length of season in days
Pstrong  = 9.; % Precipitation in strong (wet) state, mm/day
Pweak    = 0.; % Precipitation in weak (dry) state, mm/day
tau 	   = 17; % Memory length in time steps (days) 
prmax   = 0.8; % Maximum probability of either state
p_init = 0.75; % Initial probability of strong state - p_init <= prmax
% --- Random numbers
% % Model results will differ from run to run because of the stochastic 
% % nature of the model. For robust results, average over many runs. To
% % obtain reproducable results, however, a constant stream of random 
% % numbers must be used - for this purpose, uncomment the following 
% % lines: 
% s = RandStream.create('mt19937ar','seed',5489);
% RandStream.setDefaultStream(s);

% --- Compute ---
Psum = zeros(1,nyears);
% --- Loop over years (seasons)
for i = [1:nyears]
    P = zeros(1,lseason);
    n=1;
	% --- Loop over days of season
    while n <= lseason
        pr = rand; % draw random number 0 <= pr <= 1
		% --- Determine probability
        if n>tau
            p = (sum(P(n-tau:n-1))/tau - Pweak)/(Pstrong-Pweak); % ...from memory effect
        else
            p = p_init; % ...constant in the beginning of the season
		end
		% --- Limit high and low probabilities
		if p > prmax
			p = prmax;
		elseif p < (1-prmax)
			p = (1-prmax);
		end
		% --- Assign precipitation
        if pr < p
            P(n) = Pstrong; 
        else
            P(n) = Pweak; 
        end
        n=n+1;
    end
	% --- Seasonal mean
    Psum(i) = sum(P)/lseason;
end

% --- Plot ---
bins = 0:0.2:10; % Intervals for histogram bins
figure;
hist(Psum,bins)
xlim([-2,12])
ylim([0 500])
title('day-to-day model')
xlabel('mm/day')

% --- END ---