% Script to project Hispanic population in Redding.
% Input the Hispanic Data
% The data for Redding starts out in a matrix called "data".
hispanic = data(:,4); %The hispanic data is in the 4th column.
year = data(:,10); %The year is in the 10th column.
% The 2040 projection is in the 8th row.
hispanic2040 = hispanic(8); %Store the 2040 projection elsewhere.
hispanic(8)=[]; %Delete the 2040 projection from the hispanic vector.
year(8)=[];

% Plot the historical population data.
subplot(2,1,1)
plot(year,hispanic,'.')
xlabel('Year')
ylabel('Population')
title('Hispanic population of Redding')

% GOMPERTZ PROJECTION
% I will use 75000 as my ceiling
c = 75000;

% Perform input evaluation:
top = log(hispanic(3:7))-log(hispanic(2:6)); %log(P(t+1))-log(P(t))
bottom = log(hispanic(2:6))-log(hispanic(1:5)); %log(P(t))-log(P(t-1))
z = top./bottom
%Plot z value.  It should be relatively constant if the data follows a
%Gompertz trend.
subplot(2,1,2)
x = 1:5 ; %There are only 5 differences
plot(x,z,'.');
title('Growth rate of the increments')

%Perform least squares regression
y = log(log(c)-log(hispanic)); %First linearize data
OLSgomp = regres(y,year-1970); %Perform OLS regression of y on year.
% I subtracted 1970 from the year because it avoids taking the exponent of
% very large numbers.
%Calculate Gompertz coefficients
%Remember, OLSgomp.b is a 2x1 vector, with the first number being the
%intercept and the 2nd number is the slope.
b0 = exp(exp(OLSgomp.b(1)));
b1 = exp(OLSgomp.b(2));
predictyear = [1970:1:2040]; % Create a vector of years from 1970-2040 for prediction.

%Predict hispanic population
% Can do this do ways, either by
%   1. predicting yhat linear coefficients and transforming to population
%   2. predicting population from Gompertz coefficients.
%Method 1
% Predict yhat
var_e = sum(OLSgomp.resid.^2)/(5); % this line and the next five lines calc components of CI95 for yhat ;
xbar = mean(year);
num = predictyear-xbar;
denom=sum((year-xbar).^2);
se_y = sqrt(var_e*((1/6)+num/denom));
t = tinv(.975,5);

yhat = OLSgomp.b(1) + OLSgomp.b(2)*(predictyear-1970); %y = m+b*x;
yhat_l= OLSgomp.b(1) + OLSgomp.b(2)*(predictyear-1970)-t*se_y ; % CI95 lower for yhat ;
yhat_u= OLSgomp.b(1) + OLSgomp.b(2)*(predictyear-1970)+t*se_y ; % CI95 upper for yhat ;

%Transform yhat to population
prediction1 = exp(log(c)-exp(yhat));
predict_l=exp(log(c)-exp(yhat_l));
predict_u=exp(log(c)-min(exp(yhat_u),log(c)));
%Method 2;
prediction2 = c*b0.^(-b1.^(predictyear-1970));

%Create two plots of prediction, one that is linear, one that is in actual
%population.
figure; %Create new page for plots
subplot(2,1,1)
plot(predictyear,yhat)
hold on
plot(year,y,'.') %Plot actual data on top of prediction.
hold off
title('Linearized Gompertz projection of hispanic population')
xlabel('year')
ylabel('Transformed population')
subplot(2,1,2)
plot(predictyear,prediction2) %I could have used prediction1 as well
hold on
plot(year,hispanic,'.') %Plot actual data on top of prediction
plot(predictyear,predict_l,'-.r')
plot(predictyear,predict_u,'-.r')
hold off
title('Gompertz projection of hispanic population')
xlabel('year')
ylabel('Population')


