% imports weekly vaccination data and smooths it with hpfilter
% then converts it to new vaccines per 10,000 in each state
% last updated feb 24
% WV data changed for week 22622 - 22636 

clear;
clc;
clf;
x=readtable("weekly_CDC.csv");
y=readtable("state_populations.xlsx");
%note states are sorted alphabetically by abreviation, not name
state=x.Properties.VariableNames;
state=state(1,2:end);
date=x(:,1);
pop=y(:,3);
p=pop{:,:};
pt=p';
ik=6;
ccc=0;

% data from May 5, 2021 to Feb 16, 2022.
xx=x(9:end, 2:end);
xx=xx{:,:}; 

[T1,n1]=size(xx);
x1=zeros(T1,n1);
for bb=1:n1
    for cc=1:T1
        x1(cc,bb) = xx(cc,bb)./pt(1,bb)*100;
    end
end

x2=log(x1);
[Z, ~] = hpfilter(x2,400);

c=array2table(Z,'VariableNames',state);
writetable(c,'cum_vax_smoothed_0603.csv');


% now get the number of new vaccinations per 10,000 people
x3=log(xx);
[Z,~]=hpfilter(x3,400);
z3=exp(Z);
[T,n]=size(z3);

zz=zeros(T-1,n);
for jj=1:n
for ii=2:T
    zz(ii-1,jj) = z3(ii,jj) - z3(ii-1,jj);
end
end

pz=zeros(T-1,n);
for ll=1:n
    for mm=1:T-1
    pz(mm,ll)= zz(mm,ll)./pt(1,ll)*10000;
    end
end


b=array2table(pz,'VariableNames',state);
writetable(b,'newvaxper10000_smoothed_0225.csv');


