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];


tickers ={'AAPL','XOM','MSFT','BAC^I','IBM','CVX','GE','WMT','T','JNJ','PG','WFC','KO','PFE','JPM','GOOG','PM','VOD','ORCL','INTC','MRK','VZ','BRK/A','QCOM','PEP','CSCO','AMZN','ABT','MCD','SLB','C','BRK/B','BAC','RY','HD','DIS','TD','UTX','OXY','AXP','KFT','COP','MO','CAT','V','CMCSA','MMM','USB','AIG','EMC','CVS','BNS','UNH','BA','UPS','BMY','AMGN','UNP','GS','MA','EBAY','DD','HPQ','LLY','CL','SU','SPG','EPD','HON','LVS','ESRX','SBUX','NKE','ACN','F','MDT','SO','MON','TEVA','DOW','GILD','ABX','IMO','LOW','DHR','CNI','DB','TGT','PCLN','BMO','POT','MET','COST','EMR','BIDU','TXN','GM','TWX','CNQ','FCX','APC','PNC','PX','APA','BP','YUM','DTV','DE','NWSA','QQQ','MS','BCE','ENB','BIIB','TJX','KMB','CELG','BAX','TRP','NOV','HAL','D','GG','CM','MCI','WAG','RDS/A','DUK','LMT','DELL','EOG','AGN','FDX','ERIC','BK','SCCO','ITW','NEE','DVN','ADP','AMT','COV','INFY','ACE','BLK','EXC','TYC','BEN','PCP','TRV','COF','KMI','PRU','GIS','TWC','TRI','CVE','GD','PSA','AMX','CTL','VIAB','NSC','RAI','CSX','MFC','LYB','NEM','ISRG','STT','BBT','WLP','JCI','ADM','MPV','CTSH','MCK','CRM','COH','EP','GLW','CMI','SYK','AFL','CBS','BAM','TCK','ASML','CB','TMO','SE','FE','WMB','CCL','KMP','PSX','MRO','EQR','PCG','AEP','MMC','ECL','YHOO','SHPGY','CMCSK','APD','BHI','K','RTN','GSK','DFS','WPZ','SDRL','HES','ED','LO','SCHW','M','MJN','CME','HNZ','BRCM','PBR/A','VTR','WPPGY','NBL','TSM','ALL','RIG','GGP','INTU','HCP','VFC','VALE','BBBY','WFM','CCI','NVS','SYY','ALXN','L','F^A','PLD','BDX','PGN','ADBE','DG','PPG','VNO','BXP','PEG','VALE/P','NLY','PPL','AON','GR','WM','SRE','NOC','NWS','VRX','ECA','CTXS','TROW','ETN','AZO','MOS','RCI','AET','LGI','MSI','CAH','LTD','TEL','PCAR','AMAT','EL','ROST','EIX','A','GWW','CLR','IP','OMC','ATVI','AVB','EMQ','GPS','MPC','NTAP','VMW','CERN','STX','MHP','TVE','SLF','DEO','AGU','PGR','RDS/B','ORLY','PXD','XEL','ITUB','IIF','HUM','FAST','MAR','SLE','CI','KR','CMG','FITB','PH','WYNN','CA','STI','CP','PAA','VLO','IR','SHW','SWK','TDC','BAP','KSS','ESV','ABV','DLTR','CF','REGN','HCN','CAM','NUE','TLM','HOG','STJ','MMI','HST','HCA','MXF','MNST','JWN','ETR','SYMC','AMP','CHK','CHKP','HOT','RHT','NTRS','ARMH','ZMH','MAT','ADI','DOV','ROK','HSY','ETP','WU','ALTR','PAYX','MTB','WY','VIV','CPB','CAG','NLSN','SNY','IVZ','FTI','AUY','WFT','FFIV','JNPR','RL','RRC','SPLS','XRX','FMX','GPC','RYAAY','AMTD','LINTA','SWN','CXO','SLW','AA','RSG','MGA','CHRW','MUR','CBE','ROP','UN','FIS','EW','DLPH','FLR','PRGO','KGC','TOT','AGNC','DTE','LMCA','RF','BBD','NXY','BG','WPI','FISV','DAL','ABC','ICE','AES','BEAM','NE','WDC','GRMN','BPO','DGX','ETE','DO','FRX','OKE','THI','APH','BSX','MYL','XLNX','CIE','CCJ','CPN','MCO','CLX','CCE','SNDK','HIG','TU','CNP','BWA','SJM','LULU','ADSK','DPS','LH','WEC','EXPD','KORS','BCR','TIF','KLAC','AVP','CLF','SIAL','KSU','MAC','MRVL','SJR','CNA','SIRI','VRSK','AME','MXIM','FOSL','FDO','KIM','LIFE','BTU','DVA','DLR','MMP','VRTX','RAX','AEE','IVN','PFG','NTES','BHP','LINE','COL','AVGO','CIT','DISCA','MPEL','CNX','HRL','GOLD','OKS','NVDA','KEY','WYN','BF/B','WAT','EQIX','PWE','UAL','CHD','JCP','FMC','EQT','BBY','EGO','NUAN','CE','LBTYA','COG','KMR','SRCL','LLTC','SLG','VAR','MBT','TXT','LLL','EEP','S','NI','SLM','JOY','EPB','ARG','EMN','TSCO','DNR','KMX','BAC^L','HSIC','NVO','MKC','XL','PLL','WIN','SNI','AAP','LNC','BMC','NYX','JBHT','TAP','AEM','DISH','HTZ','DRI','PETM','VRSN','MCHP','MU','FRT','TDG','BX','CLB','TRMB','NU','TV','BLL','ADS','HBC','VMED','WSH','RE','UNM','GRPN','CMA','IHS','PCL','RIMM','LUV','CHTR','BUD','HFC','FLS','MGM','FNV','CFN','WLL','BVN','UDR','NXPI','LUK','AWK','BTE','MWE','SCG','SBAC','PVH','CBG','SHLD','CMS','CINF','RCL','ANSS','AKAM','VHI','Y','SXCI','ALB','NYB','LBTYK','HBAN','HSP','ST','XRAY','RYN','AMG','EFX','ILMN','HLF','AU','PII','CPT','TIBX','LVLT','XEC','ESS','OGE','ULTA','RIO','ALV','TSN','CCK','DISCK','TRW','NWL','WCRX','CHL','SUN','VMC','MTD','IRM','QEP','WRB','DOX','DHI','OII','TS','ACGL','O','PNW','HTR','TKR','BWP','DCI','PXP','MAS','ASH','SBH','JEC','HP','IPG','CTAS','MHK','SWY','AVT','AMD','MKTAY','LKQX','EA','PAY','GIB','ENR','EXPE','CEF','ELN','NBR','INFA','BPL','XYL','DKS','RMD','LNT','AZN','TRIP','MWV','LRCX','SWKS','LNKD','PWR','IDXX','PNRA','WHR','PRE','IFF','BF/A','TMK','IBN','VAL','BEAV','UL','HRS','FL','FDS','HAS','NFX','ARE','FLEX','WTW','TCO','RJF','IGT','WFC^L','GAS','HOLX','AGCO','LSI','KBR','JBL','MCRS','AXS','LEN','IEP','ATI','TNH','PBCT','CHSI','TSS','GRA','AJG','IAG','KOS','POM','VNO^A','MKL','PKX','CVH','TOL','ANF','SNPS','AN','HRB','APOL','MSCI','HDB','LRY','FNF','MDU','CIG','CPO','FRC','REG','PNR','CSC','RHI','TEG','WES','RKT','ARW','COO','RDC','CBI','ROC','SM','BIP','BMRN','OC','HUB/B','RNR','DDR','NDAQ','MELI','SIG','URBN','LECO','IT','NFLX','NVR','NGLS','LEA','WLT','NGG','RGA','RS','BRE','BTI','OIS','BOKF','NGD','RAH','FTNT','UA','SXL','ENDP','SAI','QGEN','AEO','GNC','SI','OCR','GMCR','NFG','WCN','AFG','FTR','NS','DNKN','SPW','NVE','ARBA','TE','X','DRE','WSM','SPN','LPLA','WBC','TW','PHM','BRO','ITC','TNB','YNDX','GM^B','WLK','SINA','EDU','EGX','HCS^B','IACI','PDCO','ONNN','RGP','WAB','NUC','FBHS','CREE','NCR','TAC','G','AMLN','GPRO','UHS','DNB','ZION','OI','MLM','SNA','DRC','C^P','KEX','TPX','YPF','RGLD','CFR','TEF','ATR','HCBK','WR','IEX','HUN','ALSN','CBSH','GPN','SNH','SPIL','NRG','EGN','STR','ARCC','JLL','ATML','CX','SWI','APU','WTM','MSM','CNQR','TWTC','STZ','RPM','SEE','PBI','SEIC','SOA','GFI','HPT','TM','ASNA','CSL','PENN','VR','NDSN','FLIR','ACC','TUP','ERIE','TSLA','PHI','ACAS','CLH','ABB','CPRT','C^H','AVY','FLT','LM','EWBC','WPX','HCC','AIV','KRC','GNTX','KMT','WRI','LAZ','UGI','SFD','SON','JAH','AIZ','SMG','MD','NATI','ERJ','BMS','GIL','BSBR','HAR','MAN','DEI','ERF','NVLS','ANR','SID','MEOH','SLH','GGG','VMI','TGI','GCI','WTR','BSAC','PGH','GDI','JEF','PKI','ROL','TSO','SEP','FMCN','CDNS','NETC','CFX','CRI','LPI','UFS','HME','KKR','N','SLGN','FNFG','AGP','LEG','TER','URS','CBD','BMR','RVBD','DDS','PDM','EV','SPLK','SBNY','FSL','BWC','SPR','ROVI','SD','CIM','TFSL','C^I','OTEX','GME','ATO','UPL','IM','ELS','SLXP','NNN','CTRP','CVC','HGSI','EXR','SKT','JAZZ','VIP','ETFC','GNW','INCY','NEU','CYT','JKHY','HTS','AWH','WCC','CXW','TLK','MORN','OAS','NUS'};

%top 75 by capitalization with complete records
tickers = {'XOM'    'MMM'    'CNP'};


start_date='01011982';
end_date='15122012';
end_date='01052013';
stocks = hist_stock_data(start_date,end_date,tickers{:}, 'frequency','m');
nyears=32;

n0 = length(stocks);
dta = NaN(1,(nyears*12-7));
iok=1;
for ii=1:n0
  wrk = flipud(stocks(ii).Close);
  if length(wrk)==(nyears*12-7)
    dta(iok,:) = wrk;
    sname{iok}=stocks(ii).Ticker;
    iok=iok+1;
  end
end

rdta = log(dta(:,5:12:end)./dta(:,4:12:end));

n=size(rdta,1)
for ii=1:n
  c(ii)=corr2(air(:,2),rdta(ii,1:29)');
end



[b,se,pval,inmodel0,stats,next,hist]=stepwisefit(rdta(:,1:29)',sqrt(air(:,2)));

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

inmodel = inmodel0;% & (pval'<0.01) ;

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

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')

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)]';

