clear
air =[ 1982 735.403301816046
1983 955.709541669618
1984 836.702280403177
1985 759.837968749349
1986 743.081159793134
1987 697.243120496166
1988 961.544289566035
1989 866.734443173754
1990 908.628765340288
1991 785.456208433496
1992 785.023827180553
1993 866.907336535751
1994 952.85246593786
1995 800.380205594926
1996 858.673272277702
1997 871.561407987435
1998 859.468119732363
1999 822.134574790892
2000 778.44088008849
2001 793.614203610427
2002 662.197796567331
2003 849.679313577223
2004 744.987582369206
2005 857.190495548418
2006 870.222630405468
2007 935.418611378442
2008 887.323177948036
2009 667.971028328757
2010 845.966928408054];

nyears=32;

for im=1:10000
rdta = randn(75,32);

[b,se,pval,inmodel0,stats,next,hist0]=stepwisefit(rdta(:,1:29)',sqrt(air(:,2)),'penter',0.01,'display','off');

nnz(inmodel0);
p=pval(inmodel0);
find(inmodel0);

inmodel = inmodel0 & (pval'<0.01) ;

if nnz(inmodel)

y = sqrt(air(:,2));
x = [rdta([inmodel],1:29);ones(1,29)]';
hat = x*inv(x'*x);
b0=regress(y,x);
yhat = x*b0;

y = sqrt(air(:,2))';
x = rdta([inmodel],1:29);
xm = mean(x,2);
xa = x - xm(:,ones(1,29));
ym = mean(y);
ya = y - ym;

yhat_cv = NaN(29,1);
for ii=1:29
  ind = setdiff([1:29],ii);
  b = ya(ind)*xa(:,ind)'*inv(xa(:,ind)*xa(:,ind)');
  yhat_cv(ii) = b*xa(:,ii)+ym;
end

if 0
subplot(3,6,[1:3,7:9,13:15])
h=plot(y.^2,yhat_cv.^2,'bo');
set(h,'MarkerFaceColor','k')
hold on
plot([650:1000],[[650:1000]',ones(351,1)]*regress(y'.^2,[yhat_cv'.^2,ones(29,1)]),'--k');
hold off
grid
xlabel('AIR')
ylabel('AIR prediction CV')
text(665,930,['Cross-validated correlation = ',num2str(corr2(y'.^2,yhat_cv.^2),2)])
text(665,910,['Cross-validated RMSE = ',num2str(sqrt(mean((y'.^2-yhat_cv.^2).^2)),3),'mm'])
text(665,890,['Climatology std. = ',num2str(sqrt(mean((y'.^2-ym.^2).^2)),3),'mm'])

subplot(3,6,[4:6,10:12])
bar([1982:2010],[y.^2-mean(y.^2);yhat_cv.^2-mean(yhat_cv.^2)]')
set(gca,'xlim',[1981,(nyears-29+2010+1)],'ylim',[-175,175])
grid

xout = [rdta([inmodel],[30:nyears]);ones(1,nyears-29)]';
yout = xout*b0;
hold on
h=plot([2011:(nyears-29+2010)],[yout.^2-mean(yhat.^2)],'bs')
set(h,'MarkerFaceColor','r')
plot([(nyears-29+2010)],[yout(end).^2-mean(yhat.^2)],'bx')
hold off
legend('AIR','AIR prediction CV')
end
C(im) = corr2(y'.^2,yhat_cv.^2);
else
  C(im) = 0;
end
end
clf
hist(C,20)
mean(C>0.59)
return


[u,s,v] = svd(rdta(abs(c)>0.2,:),'econ');
[b,se,pval,inmodel,stats,next,hist]=stepwisefit(v(1:29,:),air(:,2));
x = [v(1:29,inmodel)';ones(1,29)]';

