load moscow_smooth.dat % file with smoothed 1881-2009 Moscow July temperature data from GISS
sernumber = 100000; % number of Monte Carlo realisations
trendline = (moscow_smooth(:,2))/1.55; % create trendline normalised by std deviation of interannual variability
excount=0; % initialise extreme counter
for i=1:sernumber % loop through individual realisations of Monte Carlo series
    t = trendline + randn(129,1); % make a Monte Carlo series of trendline + noise
    tmax=-99;
    for j = 1:129
      if t(j) > tmax; tmax=t(j); if j >= 120; excount=excount+1; end; end % count heat records in last decade
    end
end
expected_records = excount/sernumber % expected number of records in last decade (ca. 0.41)
probability_due_to_trend = 100*(expected_records-0.079)/expected_records % probability for record happening due to trend (ca. 81%)