%% Modeling World Population Growth
%% Introduction
% In 1968 Paul Erlich published The Population Bomb, in which he predicted that world population would grow quickly during the 1970s, that agricultural production could not keep up, and that mass starvation in the next two decades was inevitable (see https://en.wikipedia.org/wiki/The_Population_Bomb).  As someone who grew up during those decades, I am happy to report that those predictions were wrong.
% But world population growth is still a topic of concern, and it is an open question how many people the earth can sustain while maintaining and improving our quality of life.
% In this case study, we use simple models to explain world population growth since 1950 and generate predictions for the next 50--100 years.
% For background on world population growth, watch this video from the American Museum of Natural History: https://www.youtube.com/watch?v=PUwmA3Q0_OE.

clear all

%
%% Data
% The Wikipedia article on world population estimates (https://en.wikipedia.org/wiki/World_population_estimates) summarizes data from the United Nation and U.S. Census Bureau, which I have converted to CSV files.
% Here's the data from the UN Department of Economic and Social Affairs (DESA):

filename = '../data/world_pop_un.csv';
m = readmatrix(filename);
year = m(:,1);
un = m(:,2)'

% And from the U.S. Census Bureau:

filename = '../data/world_pop_census.csv';
m = readmatrix(filename);
census = m(:,2)'

% Here's what their estimates look like from 1950 to 2016.

clf; hold on

plot(year, un, ':', 'LineWidth', 2)
plot(year, census, '--', 'LineWidth', 2)
legend('UN DESA', 'US Census', 'Location', 'northwest')

xlabel('Year')
ylabel('World population (billion)')

% To get a sense of how accurate these estimates are, we can see how different they are.
% The following expression computes the elementwise differences between the two series, then divides through by the UN value to produce relative errors (see https://en.wikipedia.org/wiki/Approximation_error), then finds the largest element.
% The largest relative error between the estimates is about 1.3%.

rel_error = abs(census - un) ./ un
percent_error = rel_error * 100
max_percent_error = max(percent_error)

%
%% Constant growth model
% Suppose we want to predict world population growth over the next 50 or 100 years.  We can do that by developing a model that describes how populations grow, fitting the model to the data we have so far, and then using the model to generate predictions.
% I will demonstrate this process starting with simple models and gradually improving them.
% Although there is some curvature in the plotted estimates, it looks like world population growth has been close to linear since 1960 or so.  So we'll start with a model that has constant growth.
% To fit the model to the data, I'll compute the average annual growth from 1950 to 2016.  Since the UN and Census data are so close, I'll use the Census data.
% Here are the first and last elements of the series (in billions):

p_0 = census(1)
p_end = census(end)

% And here is the elapsed time in years:

t_0 = year(1)
t_end = year(end)
elapsed_time = t_end - t_0

%  Then we can compute the average annual growth in billions of people per year.

total_growth = p_end - p_0
annual_growth = total_growth / elapsed_time

% Our simulation will start with the observed population in 1950, p_0, and add annual_growth each year.  To store the results, I'll create a vector the same size as census.

results = zeros(size(census));

%  Initially results is all zeros.  I'll initialize it so the starting value, in 1950, is the 1950 population estimated by the US Census.

results(1) = census(1)

% Now I'll run a loop from 1 to the length of census, filling in the elements of results.

for i = 1:length(census)-1
    results(i+1) = results(i) + annual_growth;
end

% Here's what the results looks like, compared to the actual data.

clf; hold on

plot(year, un, ':', 'LineWidth', 2)
plot(year, census, '--', 'LineWidth', 2)
plot(year, results, 'LineWidth', 2)
legend('UN DESA', 'US Census', 'model', 'Location', 'northwest')

xlabel('Year') 
ylabel('World population (billion)')
title('Constant growth model')

% The model does not fit the data particularly well from 1950 to 1990, but after that, it's pretty good.  Nevertheless, there are problems:
% There is no obvious mechanism that could cause population growth to be constant from year to year.  Changes in population are determined by the fraction of people who die and the fraction of people who give birth, so we expect them to depend on the current population.
% According to this model, we would expect the population to keep growing at the same rate forever, and that does not seem reasonable.
% We'll try out some different models in the next notebook.