%  This script calculates extreme values within a week from an AR(1)
%  red-noise process.

clear

%  Start by generating a red noise process
ndays = 10000;  %  Length of record
rho1 = 0.7;     %  lag-1 autocorrelation - change this

%  Initialize array
y = repmat(NaN, [ndays 1]);
y(1) = randn(1);

%  Now, generate red noise process
for i = 2:ndays;
  y(i) = rho1*y(i-1) + sqrt(1-rho1^2)*randn(1);
end

mean(y)    %  This should be close to zero
var(y)     %  This should be close to one

%  OK, let's go through and identify extreme values
days_per_week = 7;   %  Change this around if you'd like.
nweeks = floor(ndays / days_per_week);

%  Now, initialize our counting array:  the first row will be the number
%  of times the maximum value of 'y' occurs on a given day of the week,
%  and the second row will be for the mininum value of 'y'.
day_extreme_occurs = zeros(2,days_per_week);

%  And finally, count:
for i = 1:nweeks
  ind = days_per_week*(i-1)+[1:days_per_week];
  ind1 = find(y(ind) == max(y(ind)));
  day_extreme_occurs(1,ind1) = day_extreme_occurs(1,ind1)+1;
  ind2 = find(y(ind) == min(y(ind)));
  day_extreme_occurs(2,ind2) = day_extreme_occurs(2,ind2)+1;
end

%  Finally, plot this up as a histogram
figure_tall(1); clf;
subplot(2,1,1);
  bb = bar(1:days_per_week, day_extreme_occurs'/nweeks);
  %  The 'max' bar is on the left, the 'min' bar is on the right.
  t1 = title(['\bfOccurrence of extreme values of AR(1) process ' ...
              '(relative to week):  rho = 0.7']);
  set(t1, 'fontsize', 14);
  xl = xlabel('Day of the Week', 'fontsize', 12);
  yl = ylabel('Frequency', 'fontsize', 12);
  set(gca, 'XLim', [0.5 days_per_week+0.5]);
  grid on;
  

  
%  Here's a hint for the proof

%  What is the variance of departures from the weekly mean,
%  as a function of the day of the week?
yy = y(1:(nweeks*days_per_week));
yy = reshape(yy, days_per_week, nweeks)';

%  Now, remove each week's mean prior to computing the variance
yy = yy - mean(yy, 2)*ones(1,days_per_week);

%  And, compute the variance
varyy = var(yy);

%  Now, plot up the variance as above
subplot(2,1,2);
  bb = bar(1:days_per_week, varyy);
  t1 = title(['\bfDaily variance of AR(1) process ' ...
              '(relative to week):  rho = 0.7']);
  set(t1, 'fontsize', 14);
  xl = xlabel('Day of the Week', 'fontsize', 12);
  yl = ylabel('Variance', 'fontsize', 12);
  set(gca, 'XLim', [0.5 days_per_week+0.5]);
  grid on;