%% 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);
years = m(:,1)'
un = m(:,2)'

% And from the U.S. Census Bureau:

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

%
%% 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):

system.p_0 = census(1);
system.p_end = census(end);

% And here is the elapsed time in years:

system.t_0 = years(1);
system.t_end = years(end);
system.elapsed_time = system.t_end - system.t_0;

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

total_growth = system.p_end - system.p_0
system.annual_growth = total_growth / system.elapsed_time;
system.years = years;
system

%

results = run_simulation(system, @constant_growth_update);

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

data.year = years;
data.un = un;
data.census = census;
data
results

plot_results(results, data)

%
%
%

function res = constant_growth_update(t, y, system)
    % Compute the population next year.
    
    % t: current year
    % pop: current population
    % system: structure containing parameters of the model
    
    % returns: population next year
    res = y + system.annual_growth;
end

function res = update_func2(t, y, system)
    % Compute the population next year.
    
    % pop: current population
    % t: current year
    % system: system object containing parameters of the model
    
    % returns: population next year
    births = system.birth_rate * y
    deaths = system.death_rate * y
    res = y + births - deaths
end


function res = run_simulation(system, update_func)
    results = zeros(size(system.years))
    results(1) = system.p_0;
    for i = 1:length(results)-1
        t = system.years(i);
        results(i+1) = update_func(t, results(i), system);
    end
    res = results
end


function res = plot_results(results, data)
    clf; hold on

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

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