Malthusian population growth model

Data and packages

library(tidyverse)
library(readxl)
library(forecast)

setwd("C:/Users/DELL/Desktop/Hicheel/mongoliin ediin zasag/biy daalt/population increase rate model")

popul<-read_xlsx("rus pop.xlsx",sheet=1,range=cell_cols('C'),
                 col_types = c('numeric'))
names(popul)<-c('pop')

glimpse(popul)

## Observations: 59
## Variables: 1
## $ pop <dbl> 119.897, 121.236, 122.591, 123.960, 125.345, 126.745, 127.468, ...

attach(popul)
rus_pop<-ts(pop,start = c(1960,01),frequency = 1)

Malthusian descrite model

Equation: \(p_t=(1+a)*p_{t-1}\)

p_lag<-stats::lag(rus_pop,1)
x<-ts(p_lag[-length(p_lag)],start = c(1961,1))
y=ts(rus_pop[-1],start = c(1961,1))

model1<-tslm(y~x-1)
summary(model1)

## 
## Call:
## tslm(formula = y ~ x - 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0859 -0.4597  0.2561  0.5346  1.0429 
## 
## Coefficients:
##    Estimate Std. Error t value Pr(>|t|)    
## x 1.0028486  0.0005929    1691   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6317 on 57 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 2.861e+06 on 1 and 57 DF,  p-value: < 2.2e-16

s_id<-length(rus_pop)
e_id<-s_id+6

for(i in s_id:e_id){
  a<-predict(model1,newdata=data.frame(x=tail(pop,1)))
  pop[i]<-a}

fcast<-forecast(model1,newdata = data.frame(x=tail(pop,7)))

autoplot(rus_pop)+autolayer(fcast,series='Descrite model',PI=T)+
  ylab('Population, million people')+ggtitle('Malthusian discrete model')

Malthusian continuous model

Equation: \(P_t=c*e^{at}\)

y=ts(log(popul$pop),start=c(1960,1))
x=ts(1:length(y),start = c(1960,1))

model2<-tslm(y~x)
summary(model2)

## 
## Call:
## tslm(formula = y ~ x)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.070843 -0.031036 -0.002782  0.032019  0.058356 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.8546978  0.0096144 504.938  < 2e-16 ***
## x           0.0027784  0.0002787   9.969 4.22e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03646 on 57 degrees of freedom
## Multiple R-squared:  0.6355, Adjusted R-squared:  0.6291 
## F-statistic: 99.38 on 1 and 57 DF,  p-value: 4.217e-14

sp<-length(y)+1
ep<-length(y)+7


fcast2<-forecast(model2,newdata = data.frame(x=sp:ep))

autoplot(rus_pop)+ylab('population, million people')+ggtitle('Malthusian descrite and continuous model')+
  autolayer(exp(fcast2$mean),series='Continuous model')+
  autolayer(fcast,series = 'Descrite model')

Model performance

Accuracy = But without test set, it is just training set accuracy

descrite<-accuracy(fcast)
continuous<-accuracy(fcast2)

table1<-rbind(descrite,continuous)
rownames(table1)<-c('Descrite','Continuous')
knitr::kable(table1)

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Descrite	0.0259795	0.6261950	0.5575373	0.0360529	0.4014868	0.9014939	0.9479580
Continuous	0.0000000	0.0358333	0.0313581	-0.0053236	0.6346598	6.7632335	0.9452948

Cross validation

descrite<-CV(model1)
continuous<-CV(model2)

table2<-rbind(descrite,continuous)
rownames(table1)<-c('Descrite','Continuous')
knitr::kable(table2)

	CV	AIC	AICc	BIC	AdjR2
descrite	0.4057989	-50.29884	-50.08066	-46.17796	0.9999797
continuous	0.0013832	-386.80759	-386.37123	-380.57498	0.6290994

Source: Forecasting:Principles and practice