DataMortalityRate=1/100*[0.002 0.006 0.03 0.08 0.15 0.6 2.2 5.1 9.3]; %Source: Imperial College London COVID-19 Response Team, March 16, 2020

Age=0:1:110;

AgeMiddle=[4.5,14.5,24.5,34.5,44.5,54.5,64.5,74.5,84.5];

v=spline(AgeMiddle,DataMortalityRate,Age);
v(v<0)=0;


%Absolute Numbers of Deaths and Exposure to Risk 2017 in Germany can be obtained from the Human Mortality
%Database www.mortality.org


ScenarioParameter=[1 0.5 0.1];
Overlap=0:0.0001:1;
BaseDeathRates=Deaths./ExposureToRisk;
CoronaDeaths=zeros(length(Age),length(ScenarioParameter));
OverallDeaths=zeros(length(Overlap),length(Age),length(ScenarioParameter));
DeathRates=zeros(length(Overlap),length(Age),length(ScenarioParameter));
MultiplierIncreaseDeathRates=zeros(length(Overlap),length(Age),length(ScenarioParameter));


for k=1:length(ScenarioParameter)

    CoronaDeaths(:,k)=ScenarioParameter(k)*ExposureToRisk.*v';

for i=1:length(Overlap)
    OverallDeaths(i,:,k)=Deaths+(1-Overlap(i))*CoronaDeaths(:,k);
    DeathRates(i,:,k)=OverallDeaths(i,:,k)'./ExposureToRisk;
    MultiplierIncreaseDeathRates(i,:,k)=(DeathRates(i,:,k)'./BaseDeathRates);
end


end

%Additional Deaths:
NewDeathsAge=zeros(length(Overlap),length(Age),length(ScenarioParameter));
AgeBins=[10 20 30 40 50 60 70 80 90];
NewDeathsBins=zeros(length(Overlap),length(AgeBins),length(ScenarioParameter));
DeathBins=zeros(length(AgeBins),1);
for k=1:length(ScenarioParameter)
    for i=1:length(Overlap)
        NewDeathsAge(i,:,k)=(1-Overlap(i))*CoronaDeaths(:,k);
    end
    for j=1:length(AgeBins)
        if(j<length(AgeBins))
            NewDeathsBins(:,j,k)=sum(NewDeathsAge(:,AgeBins(j)-9:AgeBins(j)+1,k),2);
            DeathBins(j)=sum(Deaths(AgeBins(j)-9:AgeBins(j)+1));
        else
            NewDeathsBins(:,j,k)=sum(NewDeathsAge(:,AgeBins(j)-9:length(Age),k),2);
            DeathBins(j)=sum(Deaths(AgeBins(j)-9:length(Age)));
        end
    end
end


%Age Increase Effect:
NewAge=zeros(length(Overlap),98-26+1,length(ScenarioParameter));
AgeDiff=zeros(length(Overlap),98-26+1,length(ScenarioParameter));
for k=1:length(ScenarioParameter)
    for i=1:length(Overlap)
        for j=26:1:98
            [minValue,closestIndex]=min(abs(BaseDeathRates(j:end)-DeathRates(i,j,k)));
            NewAge(i,j-26+1,k)=Age(closestIndex+j-1);
            AgeDiff(i,j-26+1,k)=NewAge(i,j-26+1,k)-Age(j);
        end
    end
end


%Period Life Expectancy Calculation from: 
%Statistisches Bundesamt Sterbetafel 2016/2018
%- Methoden- und Ergebnisbericht zur laufenden Berechnung von
%Periodensterbetafeln für Deutschland und die Bundesländer 5. November 2019

lx=zeros(1,length(Age));
lx(1)=100000;
for i=1:length(Age)-1
    lx(i+1)=lx(i)*(1-BaseDeathRates(i));
end

Lx=zeros(1,length(Age));
for i=1:length(Age)-1
    Lx(i)=1/2*(lx(i)+lx(i+1));
end

Tx=zeros(1,length(Age));
for i=1:length(Age)
    Tx(i)=sum(Lx(i:end));
end

exRegular=Tx./lx;

%Calculate Decreased Period Life Expectancy due to COVID: 
exCorona=zeros(length(Overlap),length(Age),length(ScenarioParameter));
for k=1:length(ScenarioParameter)
    for o=1:length(Overlap)
        for j=1:length(Age)
            lx=zeros(1,length(Age));
            lx(1)=100000;
            for i=1:length(Age)-1
                if(i==j)
                    lx(i+1)=lx(i)*(1-DeathRates(o,i,k));
                else
                    lx(i+1)=lx(i)*(1-BaseDeathRates(i));
                end
            end

            Lx=zeros(1,length(Age));
            for i=1:length(Age)-1
                Lx(i)=1/2*(lx(i)+lx(i+1));
            end

            Tx=zeros(1,length(Age));
            for i=1:length(Age)
                Tx(i)=sum(Lx(i:end));
            end
            exh=Tx./lx;
            exCorona(o,j,k)=exh(j);
        end
    end
 end

ExpectancyDiff=zeros(length(Overlap),length(Age),length(ScenarioParameter));
for k=1:length(ScenarioParameter)
    for i=1:length(Age)
        ExpectancyDiff(:,i,k)=exRegular(i)-exCorona(:,i,k);
    end
end


