%BetweenYrCorr: this program generates sets of within-year, 
% auto- and cross-correlated (with one time step) 
% vital rates. 
clear all;

%*****************Simulation Parameters***********************
% parameters for 3 vital rates: 
% a beta, a lognormal, and then a stretched beta.
vrmeans= [0.77,2,5]; % means
vrvars=	[0.02,0.5,2]; % variances (not standard deviations)
% minimum and maximum values for each vital rate: 
% zeros are place holders for rates that are not stretched betas 
vrmins= [0 0 0];
vrmaxs= [0 0 24];
% positive semi-definite matrix for within year correlations
% (corrected if original was not good):
corrin=...
[1 		0.3 	-0.2;
0.3		1		0.6;
-0.2		0.6		1];

% then the auto- and cross-correlations for one step. 
% The form should be columns of v1(t), v2(t), etc...
% and rows of v1(t+1), v2(t+1), etc.. where v1(t) is
% vital rate 1 in year t
corrout=[...
0.5	0.3	-0.1;
0.3	0.2	0.3;
-0.1	0.3	0.7];

yrspan = 50; 
 % this is the number of years of correlation info to use to 
 % simulate the correlation pattern --
 % more years are better. 50 is quite accurate. 
tmax = 50; 	% number of years of vital rates to simulate;
%*************************************************************

np = length(vrmeans);
results = [];
randn('state',sum(100*clock)); % seeds random numbers

%--------creating and using the big correlation matrix, M----
% this set of loops makes the big correlation matrix (M)
% with multi-year correlations: the  if statements are used to 
% estimate the correct correlations with increasing time lags, 
% always assuming that all long-time-lag correlations are only 
% caused by within-year and one-time-step correlations
for ii = 1:yrspan;   
 for jj=1:yrspan;
  		if ii==jj, litmx = corrin;
  			else litmx =corrout; end;
		expo=1;
		if ii > jj, expo = ii-jj; litmx = (litmx)^expo; end;
		if ii < jj expo = jj-ii;  litmx = (litmx')^expo; end;
 		for ip=1:np;
 		  for jp = 1:np;
		  M(ip+np*(ii-1),jp+np*(jj-1)) = litmx(ip,jp);
		  end; % jp
		end; % ip
 end; % jj
end; % ii

% get the eigenvalues for calculating the M12 matrix 
% (here called zfull):
[W,D] = eig(M); % getting the eigenvalues and vectors

% now, a check for negative eigenvalues -- if you have them,
% it sets negatives and small positive values = 0 
check = min(min(D));	% are the smallest eigenvalues negative?
if check < 0
	disp('The min eigenvalue is < 0. Eigenvalues are:')
	disp(diag(D))
	disp('hit enter to continue with approximation')
	pause
	maxneg = max(max(abs(D(find(D<=0))))); 
	% maxneg is the largest negative eigenvalue
	D(find(abs(D<=maxneg))) = 0; % set the negatives = 0
	disp('Corrected eigenvalues are:')
	disp(diag(D))
	newfullmx = W*newD*W'; % make a corrected matrix
	for ii=1:np % change from covariances to correlations  
	 for jj=1:np
	if newfullmx(ii,ii) ==0 | newfullmx(jj,jj)==0; 
			newfullmx(ii,jj) =0;
		else
				newfullmx(ii,jj) = ...
newfullmx(ii,jj)/((newfullmx(ii,ii)*newfullmx(jj,jj))^0.5);
		end;
	 end;
	end;
	[W,D] = eig(newfullmx);  
end; %check < 0

M12 = W*abs(D.^0.5)*W';	% the M^(1/2) matrix  
sz = length(M12); % the total number of lines in M12
	
% get the lines from the middle of M12 
% to use to generate correlations
startcase = (round(yrspan/2)*np +1);  			
zvalold=real(M12(startcase:(startcase+np-1),:));
zvalnew=real(M12((startcase+np):(startcase+2*np-1),:));
clear M12 W D; % clearing memory  
% zvalold and zvalnew are each one year of rows in M12 
% to start the whole thing off, calculate a first set of
% normals, then multiply with zvalold to get correlated normals:
newns = randn(sz,1);	
oldxy = zvalold*newns;
%-----end of: creating and using the big correlation matrix------

normresults = [];  vitalresults = [];

for tt=1:tmax % a loop to make multiple sets of rates
	disp('time is'); disp(tt);
	%update the uncorrelated random normals
	newns = [newns((np+1):sz); randn(np,1)]; 
	% make the new set of correlated normals	
	newxy= zvalnew*newns;						
	normresults = [normresults; oldxy', newxy']; % save results 
	oldxy = newxy;	% make the new correlated rates old							 

% now convert correlated normals to the correct distributions:
	vrate1 = ...
	 betaval(vrmeans(1),(vrvars(1))^0.5,stnormfx(newxy(1)));
    vrate2 = lnorms(vrmeans(2),vrvars(2),newxy(2));
    vrate3 = stretchbetaval(vrmeans(3),(vrvars(3))^0.5,...
				vrmins(3), vrmaxs(3),stnormfx(newxy(3)));
	vitalresults = [vitalresults;vrate1,vrate2,vrate3];
end; % tt

disp('the input correlations (all rates and one step')
disp('correlations) are:');
disp([corrin corrout';corrout,corrin])
disp('the correlations of the normals are:')
disp(corrcoef(normresults));
disp('the within-year correlations of the vital rates are:');
disp(corrcoef(vitalresults));
disp('the input means and variances were');
vrmeans, vrvars
disp('means and variances of the simulated rates are:')
meanrates = mean(vitalresults)
variances = var(vitalresults)