% Ted Crone style indicator
vlog_emp=log(vemp);
CycleIndicator=vlog_emp-hpfilter(vlog_emp,1600/4^4);
CycleCorr=corr(CycleIndicator).*(corr(CycleIndicator)>0);
CycleCorr=CycleCorr-diag(diag(CycleCorr));
% prctile(reshape(CycleCorr(tril(ones(size(CycleCorr))==1,-1)),1,[]), [5 25 50 75 95])
%%
DistanceCorr=xlsread('../Data/Euclidean_Distance_Matrix.xlsx');
if size(DistanceCorr,1)>size(DistanceCorr,2)
    DistanceCorr(end,:)=[];
elseif size(DistanceCorr,1)<size(DistanceCorr,2)
    DistanceCorr(:,end)=[];
end;
DistanceMSA=DistanceCorr(2:end,1);
MySetDiff=setdiff(DistanceMSA,vMSA);
[rr,~]=find(repmat(DistanceMSA,1,size(MySetDiff,1))==repmat(MySetDiff',size(DistanceMSA,1),1));
DistanceCorr(rr+1,:)=[]; % add one to account for the fact that 1st row and column are empty
DistanceCorr(:,rr+1)=[];
if ~isequal( DistanceCorr(2:end,1), vMSA(:))
    error('Distance matrix wrong');
end;
DistanceCorr(1,:)=[];
DistanceCorr(:,1)=[];
Original_DistanceCorr=DistanceCorr;
DistanceCorr( DistanceCorr>0 )=1./(DistanceCorr(DistanceCorr>0).^2);
%%
StateCorr=csvread('../Data/StateCorr.csv');
if size(StateCorr,1)>size(StateCorr,2)
    StateCorr(end,:)=[];
elseif size(StateCorr,1)<size(StateCorr,2)
    StateCorr(:,end)=[];
end;
StateMSA=StateCorr(2:end,1);
MySetDiff=setdiff(StateMSA,vMSA);
[rr,~]=find(repmat(StateMSA,1,size(MySetDiff,1))==repmat(MySetDiff',size(StateMSA,1),1));
StateCorr(rr+1,:)=[]; % add one to account for the fact that 1st row and column are empty
StateCorr(:,rr+1)=[];
if ~isequal( StateCorr(2:end,1), vMSA(:))
    error('Distance matrix wrong');
end;
StateCorr=tril(StateCorr(2:end,2:end))+tril(StateCorr(2:end,2:end))'-eye(size(StateCorr,1)-1)*2;

%% Read additional geographical data
N=length(cMSA);
[aa,bb]=xlsread('../Data/MSA Land Area 2010.xlsx');
MSA_id=aa(:,1);
MSA_area=aa(:,3)/2589988; %m^2 to miles^2

varea=zeros(N,1);
for mm=1:N
    varea(mm,1)= MSA_area(MSA_id==vMSA(mm)) ;
end;
vdensity=vpop(1:T,:)./(ones(T,1)*varea');

[aa,bb]=xlsread('../Data/MSA_HDD_CDD_71-00.xlsx');
MSA_id=aa(:,1);
HDD=zeros(N,1);
CDD=zeros(N,1);
for mm=1:N
    HDD(mm,1)=aa(MSA_id==vMSA(mm),3);
    CDD(mm,1)=aa(MSA_id==vMSA(mm),4);
end;
