% Liying Lu
% CDS 410
% Project 3: COVID-19 Modeling and Maximum Likelihood Estimation
clear;clc

% input the data from source
T = [1:17]; % Time from day 1 to 17
Y = [9;6;15;10;16;10;18;18;29;28;14;13;4;31;27;30;52];
T1 = [1:35];
T2 = [1:50];

% code for models
LSF = "Least Squares Fit";
MLE = "Maximum Likelihood Estimation"; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Question 1: Find the best fit line by first trying a least squares fit. 
A = zeros(length(T),2);
A(:,1) = 1;
A(:,2) = T;
disp("For least square fit, a1 and a2 a1 and a2 estimated using data for T=1:17 are ")
a = A\log(Y) % a1 and a2

% The least squares fit for T1 
A = zeros(length(T1),2);
A(:,1) = 1;
A(:,2) = T1;
predicted = A * a; % predicted Log(Y)

% Figure 1: Number of cases reported daily from T=1 to T=17 with least square fit to T1 = 35.
LogScaleWithPrediction(T,Y,T1,a,LSF)

% Firgure 2: plot data on linear scale with least square fit to T1 = 35.
LinearScaleWithPrediction(T,Y,T1,a,LSF)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Question 2: Compute the standard deviation of the fit
stddev = std(log(Y) - predicted(1:17));
% Figure 3: 95% confidence interval of prediction to T1 = 35.
upper = exp(a(1) + a(2)*T + 2*stddev);
lower = exp(a(1) + a(2)*T - 2*stddev);
LogScaleWithPredConfInterval(T,Y,T1,a,upper,lower,LSF)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Question 3: Cumulative cases with error bounds out to March 11th (T2=50).
% Figure 4: 95% confidence interval of cumulative prediction to T2 = 50.
upper = cumsum(exp(a(1) + a(2)*T2 + 2*stddev));
lower = cumsum(exp(a(1) + a(2)*T2 - 2*stddev));
CumSumLogScalePredConfInterval(T,Y,T2,a,upper,lower,LSF)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Question 4: Derivative of L(a1, a2). See DF.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Question 5: Second derivative or Hessian matrix of L. See HL.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Question 6: Use multivariable Newton's method to find a1 and a2 that maximize 
% L and so are the most likely values.
stop = 100;
a0 = [2.12; 0.07]; % initial guess
for i = 1: stop
  a0 = a0 - HL(a0, T)\DF(a0, T, Y);
end
disp("The maximum likihood estimation of a1 and a2 estimated using data for T=1:17 are ")
a0  % [2.113807; 0.085167]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Question 7: Repeat the plots from q1-q3 using the new a1 and a2. 
% Figure 5: Number of cases reported daily from T=1 to T=17 with maximum likelihood estimation to T1=35.
LogScaleWithPrediction(T,Y,T1,a0,MLE)

% Firgure 6: plot data on non-log scale
LinearScaleWithPrediction(T,Y,T1,a0,MLE)

% Figure 7: 
% Poisson model with variance lambda, so the confidence intervals changes to 
upper = exp(a0(1)+a0(2)*T) + 2*sqrt(exp(a0(1)+a0(2)*T));
lower = exp(a0(1)+a0(2)*T) - 2*sqrt(exp(a0(1)+a0(2)*T));
LogScaleWithPredConfInterval(T,Y,T1,a0,upper,lower,MLE)

% Figure 8: Cumulative cases
upper = cumsum(exp(a0(1)+a0(2)*T2) + 2*sqrt(exp(a0(1)+a0(2)*T2)));
lower = cumsum(exp(a0(1)+a0(2)*T2) - 2*sqrt(exp(a0(1)+a0(2)*T2)));
CumSumLogScalePredConfInterval(T,Y,T2,a0,upper,lower,MLE)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Question 8: 
T3 = [1:51]; % Jan. 22nd to Mar. 12th
Y_cumsum = [cumsum(Y);
            354;382;461;481;526;587;608;697;781;896;
            999;1124;1212;1385;1715;2055;2429;2764;3323;4288;
            5364;6780;8559;10292;12746;14905;17873;21399;25404;29256;
            33627;38170;45421;53763];
Y2 = [Y_cumsum(1); diff(Y_cumsum) ]; % daily cases from Feb 7th to Mar.12th
% Figure 9: Add latest data to least squares fit
upper = cumsum(exp(a(1) + a(2)*T3 + 2*stddev));
lower = cumsum(exp(a(1) + a(2)*T3 - 2*stddev));
CumSumLogScalePredConfInterval(T3,Y2,T3,a,upper,lower,LSF)

% Figure 10: Add latest data to maximum likelihood estimation
upper = cumsum(exp(a0(1)+a0(2)*T3) + 2*sqrt(exp(a0(1)+a0(2)*T3)));
lower = cumsum(exp(a0(1)+a0(2)*T3) - 2*sqrt(exp(a0(1)+a0(2)*T3)));       
CumSumLogScalePredConfInterval(T3,Y2,T3,a0,upper,lower,MLE)

% Comparison between the Least Square Fit and Maximum Likelihood Fit:
% The maximum likelihood estimation has a smaller confidence interval than the 
% least square fit. Both models only perform better from T=1 to T=30. As T grows
% greater than 30, the number of cumulative cases increases at the faster rate 
% than the exponential rate. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Question 9: Refit both models using currently available data and run forecast
% out to 200 days.
T4 = [1:200]; % 200 days
% Model 1: Least Square Fit
A = zeros(length(T3),2);
A(:,1) = 1;
A(:,2) = T3;
disp("For least square fit, a1 and a2 estimated using updated data for T=1:51 are ")
a = A\log(Y2) % a1 and a2

A = zeros(length(T4),2);
A(:,1) = 1;
A(:,2) = T4;
predicted = A * a; % predicted Log(Y) for T4 = 200 days
stddev = std(log(Y2) - predicted(1:length(T3)));
% Figure 11
upper = cumsum(exp(a(1) + a(2)*T4 + 2*stddev));
lower = cumsum(exp(a(1) + a(2)*T4 - 2*stddev));
CumSumLogScalePredConfInterval(T3,Y2,T4,a,upper,lower,LSF)

% Model 2: Maximum Likelihood Estimation
stop = 100;
a0 = [1.227; 0.144]; % initial guess
for i = 1: stop
  a0 = a0 - HL(a0, T3)\DF(a0, T3, Y2);
end
disp("The maximum likihood estimation of a1 and a2 estimated using updated data for T=1:51 are ")
a0
% Figure 12
upper = cumsum(exp(a0(1)+a0(2)*T4) + 2*sqrt(exp(a0(1)+a0(2)*T4)));
lower = cumsum(exp(a0(1)+a0(2)*T4) - 2*sqrt(exp(a0(1)+a0(2)*T4)));  
CumSumLogScalePredConfInterval(T3,Y2,T4,a0,upper,lower,MLE)

% The Least Squares Fit model predicts a million cumulative cases at about day 72.
% The Maximum Likelihood Estimation predicts a million cumulative cases at about day 68.
% Based on the previous check on updated data, both models predict successfully until day 30.
% So for the updated models, they may predict successfully until day 80 depending on the rate of growth of the cases.




















