Generating the Figures in Chapter 10 of Otto and Day

Clear["Global`*"]

Figure 10.2

Clear["Global`*"]

jj[0]=0.1;aa[0]=0.9;a=0;b=10;c=0.1;d=0.2;tmax=50;
For[t=0,t<tmax,t++,{jj[t+1]=a jj[t]+b aa[t],aa[t+1]=c jj[t]+d aa[t]};]
Show[ListPlot[Table[{i,jj[i]},{i,0,tmax}],Joined->True,PlotRange->{{0,tmax},{0,50}}],ListPlot[Table[{i,aa[i]},{i,0,tmax}],Joined->True,PlotRange->{{0,tmax},{0,10}},PlotStyle->{Dashing[{0.01}]}]]
Show[ListPlot[Table[{i,((jj[i])/(jj[i]+aa[i]))},{i,0,tmax}],Joined->True,PlotRange->{{0,tmax},{0,1}}],ListPlot[Table[{i,((aa[i])/(jj[i]+aa[i]))},{i,0,tmax}],Joined->True,PlotRange->{{0,tmax},{0,1}},PlotStyle->{Dashing[{0.01}]}]]
ListPlot[Table[{i,((jj[i+1]+aa[i+1])/(jj[i]+aa[i]))},{i,0,tmax-1}],Joined->True,PlotRange->{{0,tmax},{0,2}}]
Eigenvalues[{{a,b},{c,d}}]

Show[
VectorPlot[{(a-1) x1+b x2,c x1+(d-1) x2},{x1,-1,35},{x2,-1,5},Axes->True,VectorScale->{Medium,0.5,#5&},AspectRatio->1,PlotRange->{{-1,35},{-1,5}}],
Plot[{(1/(-((-a+d+(Sqrt[(a^2)+4 b c-2 a d+(d^2)]))/(2 c)))) x,(1/(-((-a+d-(Sqrt[(a^2)+4 b c-2 a d+(d^2)]))/(2 c)))) x},{x,-1,40},AspectRatio->1,PlotRange->{{-1,35},{-1,5}},PlotStyle->{{Dashing[{0}]},{Dashing[{0}]}}],
ListPlot[Table[{jj[i],aa[i]},{i,0,tmax}],Joined->True,PlotRange->{{0,35},{0,5}},PlotStyle->Thickness[0.01]]
]

1/(1/(-((-a+d+(Sqrt[(a^2)+4 b c-2 a d+(d^2)]))/(2 c))))
1/(1/(-((-a+d-(Sqrt[(a^2)+4 b c-2 a d+(d^2)]))/(2 c))))

Figure 10.3

Clear["Global`*"]

matrix={{0,0,0,b},{sic,sii,0,0},{0,smi,smm,smr},{0,sri,srm,srr}}

sic=.92;sii=0.86;smi=0.08;smm=.8;smr=1-μ;sri=.02;srm=.19;srr=0;

tmax=200;b=.3;μ=.12;
n1[1]=1;n2[1]=2;n3[1]=1.2;n4[1]=3;
For[t=1,t<tmax,t++,
{n1[t+1],n2[t+1],n3[t+1],n4[t+1]}=matrix.{n1[t],n2[t],n3[t],n4[t]}
]
Clear[b,μ]

Show[ListPlot[Table[{i,n1[i]},{i,1,tmax}],Joined->True,PlotRange->{0,15},PlotStyle->{Dashing[{0}]}],ListPlot[Table[{i,n2[i]},{i,1,tmax}],Joined->True,PlotRange->{0,15},PlotStyle->{Dashing[{0.015}]}],ListPlot[Table[{i,n3[i]},{i,1,tmax}],Joined->True,PlotRange->{0,15},PlotStyle->{Dashing[{0.03}]}],ListPlot[Table[{i,n4[i]},{i,1,tmax}],Joined->True,PlotRange->{0,15},PlotStyle->{Dashing[{0.06}]}]]

Export["fig10-3a.EPS",%];

Show[ListPlot[Table[{i,((n1[i])/(n1[i]+n2[i]+n3[i]+n4[i]))},{i,1,tmax}],Joined->True,PlotRange->{0,1},PlotStyle->{Dashing[{0}]}],ListPlot[Table[{i,((n2[i])/(n1[i]+n2[i]+n3[i]+n4[i]))},{i,1,tmax}],Joined->True,PlotRange->{0,1},PlotStyle->{Dashing[{0.015}]}],ListPlot[Table[{i,((n3[i])/(n1[i]+n2[i]+n3[i]+n4[i]))},{i,1,tmax}],Joined->True,PlotRange->{0,1},PlotStyle->{Dashing[{0.03}]}],ListPlot[Table[{i,((n4[i])/(n1[i]+n2[i]+n3[i]+n4[i]))},{i,1,tmax}],Joined->True,PlotRange->{0,1},PlotStyle->{Dashing[{0.06}]}]]

Export["fig10-3b.EPS",%];

ListPlot[Table[{i,((n1[i+1]+n2[i+1]+n3[i+1]+n4[i+1])/(n1[i]+n2[i]+n3[i]+n4[i]))},{i,1,tmax-1}],Joined->True,PlotRange->{0,1.4}]

Export["fig10-3c.EPS",%];

Figure 10.4

matrix={{0,0,0,b},{sic,sii,0,0},{0,smi,smm,smr},{0,sri,srm,srr}}

sic=.92;sii=0.86;smi=0.08;smm=.8;smr=1-μ;sri=.02;srm=.19;srr=0;

maxeigen=Max[Eigenvalues[matrix/.μ->0.12/.b->0.3]]

Show[
ListPlot[{{0.12,maxeigen},{0.3,maxeigen}},PlotStyle->PointSize[0.02]],
Plot[Max[Eigenvalues[matrix/.b->0.3]],{μ,0.12,0.3}],
Plot[Max[Eigenvalues[matrix/.μ->0.12]],{b,0.1,0.3},PlotStyle->Dashing[{0.04}]],
PlotRange->{{0,0.31},{0.99,1.01}},
AxesOrigin->{0,1}
]

Export["Fig.eps",%];

Figure 10.5

matrix={{0,0,0,b},{sic,sii,0,0},{0,smi,smm,smr},{0,sri,srm,srr}}

sic=.92;sii=0.86;smi=0.08;smm=.8;smr=1-μ;sri=.02;srm=.19;srr=0;

tmax=50;b=.3;μ=.12;
n1[1]=10;n2[1]=10;n3[1]=10;n4[1]=10;
For[t=1,t<tmax,t++,
{n1[t+1],n2[t+1],n3[t+1],n4[t+1]}=matrix.{n1[t],n2[t],n3[t],n4[t]}
]
Clear[b,μ]

plot1=ListPlot[Table[{i,n1[i]+n2[i]+n3[i]+n4[i]},{i,1,tmax}],Joined->True,PlotRange->{0,100},PlotStyle->{Dashing[{0}]}]

n50=(n1[i]+n2[i]+n3[i]+n4[i])/.i->50

tmax=50;b=.3;μ=.12;
n1[1]=20;n2[1]=10;n3[1]=10;n4[1]=10;
For[t=1,t<tmax,t++,
{n1[t+1],n2[t+1],n3[t+1],n4[t+1]}=matrix.{n1[t],n2[t],n3[t],n4[t]}
]
Clear[b,μ]

plot2=ListPlot[Table[{i,n1[i]+n2[i]+n3[i]+n4[i]},{i,1,tmax}],Joined->True,PlotRange->{0,100},PlotStyle->{Dashing[{0.005}]}]

(n1[i]+n2[i]+n3[i]+n4[i])/.i->50

standard=(%-n50)

tmax=50;b=.3;μ=.12;
n1[1]=10;n2[1]=20;n3[1]=10;n4[1]=10;
For[t=1,t<tmax,t++,
{n1[t+1],n2[t+1],n3[t+1],n4[t+1]}=matrix.{n1[t],n2[t],n3[t],n4[t]}
]
Clear[b,μ]

plot3=ListPlot[Table[{i,n1[i]+n2[i]+n3[i]+n4[i]},{i,1,tmax}],Joined->True,PlotRange->{0,100},PlotStyle->{Dashing[{0.012}]}]

(n1[i]+n2[i]+n3[i]+n4[i])/.i->50

(%-n50)/standard

tmax=50;b=.3;μ=.12;
n1[1]=10;n2[1]=10;n3[1]=20;n4[1]=10;
For[t=1,t<tmax,t++,
{n1[t+1],n2[t+1],n3[t+1],n4[t+1]}=matrix.{n1[t],n2[t],n3[t],n4[t]}
]
Clear[b,μ]

plot4=ListPlot[Table[{i,n1[i]+n2[i]+n3[i]+n4[i]},{i,1,tmax}],Joined->True,PlotRange->{0,100},PlotStyle->{Dashing[{0.02}]}]

(n1[i]+n2[i]+n3[i]+n4[i])/.i->50

(%-n50)/standard

tmax=50;b=.3;μ=.12;
n1[1]=10;n2[1]=10;n3[1]=10;n4[1]=20;
For[t=1,t<tmax,t++,
{n1[t+1],n2[t+1],n3[t+1],n4[t+1]}=matrix.{n1[t],n2[t],n3[t],n4[t]}
]
Clear[b,μ]

(n1[i]+n2[i]+n3[i]+n4[i])/.i->1

plot5=ListPlot[Table[{i,n1[i]+n2[i]+n3[i]+n4[i]},{i,1,tmax}],Joined->True,PlotRange->{{0,50},{35,65}},PlotStyle->{Dashing[{0.06,0.025}]}]

(n1[i]+n2[i]+n3[i]+n4[i])/.i->50

(%-n50)/standard

Show[plot1,plot2,plot3,plot4,plot5,PlotRange->{{0,50},{35,65}},
AxesOrigin->{0,35}]

Export["fig10-5.EPS",%]

Figure 10.6

matrix={{m1,m2,m3,m4},{p1,0,0,0},{0,p2,0,0},{0,0,p3,0}}

matrix={{2,3,4,4},{.6,0,0,0},{0,.3,0,0},{0,0,.1,0}}

Eigensystem[matrix]

0.9767917926963138+0.21292834099431543+0.0232078518357736+0.0008431700297481771

{0.9767917926963138`,0.21292834099431543`,0.0232078518357736`,0.0008431700297481771`}/%

Eigensystem[Transpose[matrix]]

{-0.37924390001046354`,-0.4756048300067573`,-0.5711594983191037`,-0.5511360426056047`}/-0.37924390001046354`

-0.37924390001046354-0.4756048300067573-0.5711594983191037-0.5511360426056047

tmax=20;
n1[1]=1;n2[1]=2;n3[1]=1.2;n4[1]=3;
For[t=1,t<tmax,t++,
{n1[t+1],n2[t+1],n3[t+1],n4[t+1]}=matrix.{n1[t],n2[t],n3[t],n4[t]}
]

Show[ListPlot[Table[{i,n1[i]},{i,1,tmax}],Joined->True,PlotRange->All,PlotStyle->{Dashing[{0}]}],ListPlot[Table[{i,n2[i]},{i,1,tmax}],Joined->True,PlotRange->All,PlotStyle->{Dashing[{0.015}]}],ListPlot[Table[{i,n3[i]},{i,1,tmax}],Joined->True,PlotRange->All,PlotStyle->{Dashing[{0.03}]}],ListPlot[Table[{i,n4[i]},{i,1,tmax}],Joined->True,PlotRange->All,PlotStyle->{Dashing[{0.06}]}],
PlotRange->{{0,20},{0,2000}},
AxesOrigin->{0,0}]

Export["fig10-6a.EPS",%];

Show[
ListPlot[Table[{i,((n1[i])/(n1[i]+n2[i]+n3[i]+n4[i]))},{i,1,tmax}],Joined->True,PlotRange->All,PlotStyle->{Dashing[{0}]}],
ListPlot[Table[{i,((n2[i])/(n1[i]+n2[i]+n3[i]+n4[i]))},{i,1,tmax}],Joined->True,PlotRange->All,PlotStyle->{Dashing[{0.015}]}],ListPlot[Table[{i,((n3[i])/(n1[i]+n2[i]+n3[i]+n4[i]))},{i,1,tmax}],Joined->True,PlotRange->All,PlotStyle->{Dashing[{0.03}]}],
ListPlot[Table[{i,((n4[i])/(n1[i]+n2[i]+n3[i]+n4[i]))},{i,1,tmax}],Joined->True,PlotRange->All,PlotStyle->{Dashing[{0.06}]}],
PlotRange->{{0,20},{0,1}},
AxesOrigin->{0,0}]

Export["fig10-6b.EPS",%];

Canadian Demographics

Obtained from the Statistics Canada catalogues and US Census Bureau IDB 
(http://www.census.gov/ftp/pub/ipc/www/idbsprd.html)

Figure 10.7 - Total Numbers

Table 001. Total Midyear Population
---------------- --------------------
Country or area/                     
Year                       Population
---------------- --------------------
Canada

1950                       14,011,422
1955                       16,050,356
1960                       18,266,765
1965                       20,071,104
1970                       21,749,986
1975                       23,209,200
1980                       24,593,300
1985                       25,941,600
1990                       27,790,600
1995                       29,619,002
2000                       31,278,097
2005                       32,805,041
---------------- --------------------
Source: U.S. Census Bureau, International 
        Data Base.

T1={{14011422,1950},{16050356,1955},{18266765,1960},
{20071104,1965},{21749986,1970},{23209200,1975},{24593300,1980},
{25941600,1985},{27790600,1990},{29619022,1995},{31278097,2000},{32805041,2005}};

Total number of individuals

Transpose[{Transpose[T1][[2]],Transpose[T1][[1]]}];
Show[ListPlot[%,PlotStyle->PointSize[0.02]],ListPlot[%,Joined->True],PlotRange->{{1940,2010},{0,3.5 (10^7)}},
AxesOrigin->{1940,0}]

BarChart[Transpose[T1][[1]],AspectRatio->0.5,ChartLabels->Transpose[T1][[2]],BarSpacing->0.2]

"U.S. Bureau of the Census, International Data Base"

"Table 008. Vital Rates"
"--------","----------","----------","----------------","---------------","-----------"
"        ","    Births","    Deaths","   Net number of","Rate of natural","           "
"Country/"," per 1,000"," per 1,000","    migrants per","       increase","Growth rate"
"Year    ","population","population","1,000 population","      (percent)","  (percent)"
"--------","----------","----------","----------------","---------------","-----------"
"Canada"
1995    ,     13.61,      7.15,            4.54,          0.646,            1.101
"--------","----------","----------","----------------","---------------","-----------"
"Source: U.S. Bureau of the Census, International" 
"        Data Base."

Figure 10.8 - Female age distribution - 1991, 1996, & 1951

From Statistics Canada (June 1st 1991, post-censal female population estimates):

T2={{1000*(195.9+749.5),"0"},{1000*(909.7),"5"},{1000*901.5,"10"},{1000*895.9,"15"},
{1000*976.8,"20"},{1000*1168.8,"25"},{1000*1221.8,"30"},{1000*1133,"35"},{1000*1031.5,"40"},
{1000*807.5,"45"},{1000*656.6,"50"},{1000*610.1,"55"},{1000*599.4,"60"},
{1000*571.8,"65"},{1000*461.7,"70"},{1000*362.6,"75"},{1000*237.9,"80"},{1000*197.1,"85+"}};

Sum[T2[[i,1]],{i,1,18}]

BarChart[Transpose[T2][[1]],AspectRatio->0.5,ChartLabels->Transpose[T2][[2]],BarSpacing->0.2]

From Statistics Canada (July 1st 1996, female population estimates):

T2c={{954956,"0"},{984523,"5"},{987683,"10"},{976548,"15"},
{1002856,"20"},{1102079,"25"},{1297200,"30"},{1322502,"35"},{1195712,"40"},
{1074722,"45"},{833969,"50"},{670657,"55"},{616911,"60"},
{593058,"65"},{547088,"70"},{415117,"75"},{292734,"80"},{162328+87975,"85+"}};

BarChart[Transpose[T2c][[1]],AspectRatio->0.5,ChartLabels->Transpose[T2c][[2]],BarSpacing->0.2]

"Canada/1951/Total/Enu"

"Total, all ages",       14009429,        7088873,        6920556
"  0- 4         ",        1722109,         879063,         843046
"  5- 9         ",        1397825,         713873,         683952
" 10-14         ",        1130783,         575122,         555661
" 15-19         ",        1057972,         532180,         525792
" 20-24         ",        1088641,         537535,         551106
" 25-29         ",        1131215,         552812,         578403
" 30-34         ",        1042734,         512557,         530177
" 35-39         ",         999133,         503571,         495562
" 40-44         ",         868567,         445800,         422767
" 45-49         ",         744679,         387708,         356971
" 50-54         ",         662656,         340461,         322195
" 55-59         ",         570690,         292564,         278126
" 60-64         ",         506152,         264324,         241828
" 65-69         ",         433497,         228076,         205421
" 70-74         ",         315072,         160398,         154674
" 75-79         ",         188391,          94130,          94261
" 80-84         ",          96791,          45963,          50828
" 85-89         ",          39599,          17539,          22060
" 90-94         ",          10742,           4370,           6372
" 95+           ",           2181,            827,           1354

T2b={{843046,"0"},{683952,"5"},{555661,"10"},{525792,"15"},
{551106,"20"},{578403,"25"},{530177,"30"},{495562,"35"},{422767,"40"},
{356971,"45"},{322195,"50"},{278126,"55"},{241828,"60"},
{205421,"65"},{154674,"70"},{94261,"75"},{50828,"80"},{22060+6372+1354,"85+"}};

BarChart[Transpose[T2b][[1]],AspectRatio->0.5,ChartLabels->Transpose[T2b][[2]],BarSpacing->0.2]

Figure 10.9 - Mortality distribution -1991

From Statistics Canada (1991,deaths in a single year for the female population of age x):

The total number of deaths per age group is:

T3={{1344,"0"},{155,"5"},{152,"10"},{345,"15"},
{361,"20"},{556,"25"},{681,"30"},{939,"35"},{1337,"40"},
{1659,"45"},{2142,"50"},{3275,"55"},{5032,"60"},
{7763,"65"},{9757,"70"},{12949,"75"},{14261,"80"},{13381+9267+3861+906,"85+"}};

MatrixForm[Transpose[T3][[1]]]

Sum[T3[[i,1]],{i,1,18}]

BarChart[Transpose[T3][[1]],AspectRatio->0.5,ChartLabels->Transpose[T3][[2]],
BarSpacing->0.2,PlotRange->{Automatic,{-1500,30000}}]

The death rate per year individual is the number of deaths divided by 
the number of individuals:

T3a=Transpose[{N[Transpose[T3][[1]]/Transpose[T2][[1]]],
	Transpose[T2][[2]]}]

BarChart[Transpose[T3a][[1]],AspectRatio->0.5,ChartLabels->Transpose[T3a][[2]],
BarSpacing->0.2,PlotRange->{Automatic,{-0.007,0.15}}]

The death rate per census point per individual is the number of deaths divided by 
the number of individuals:

T3b=Transpose[{1-(1-N[Transpose[T3][[1]]/Transpose[T2][[1]]])^5,Transpose[T2][[2]]}]

dx=N[Transpose[T3b][[1]],10]

lx=Table[Product[1-dx[[i]],{i,1,j}],{j,0,18}]

BarChart[Transpose[T3b][[1]],AspectRatio->0.6,ChartLabels->Transpose[T3b][[2]],
BarSpacing->0.2,PlotRange->{Automatic,{-0.05,1.01}}]

Figure 10.10 - Fertility distribution - 1991

From Statistics Canada (1991, live births by age of mother):

T4={{0,"0"},{0,"5"},{265,"10"},{24180,"15"},
{80723,"20"},{150024,"25"},{107560,"30"},{33107,"35"},{4124,"40"},
{138,"45"},{0,"50"},{0,"55"},{0,"60"},
{0,"65"},{0,"70"},{0,"75"},{0,"80"},{0,"85+"}};

MatrixForm[Transpose[T4][[1]]]

The female birth rate per census point per female is the number of births divided by 
the number of individuals times 5 times the birth sex ratio:

T4c=Transpose[{N[5*0.4896418833302926*Transpose[T4][[1]]/Transpose[T2][[1]]],Transpose[T2][[2]]}]

BarChart[Transpose[T4c][[1]],AspectRatio->0.6,ChartLabels->Transpose[T4c][[2]],
BarSpacing->0.2,PlotRange->{Automatic,{-0.025,0.5}}]

Figure 10.11 - Projection - 1991 vs 1996

mx=N[5*0.4896418833302926*Transpose[T4][[1]]/Transpose[T2][[1]]]

dx=N[Transpose[T3b][[1]],10]

px=1-dx;

leslie=Join[{mx},Table[DiagonalMatrix[px][[i]],{i,1,16}],{DiagonalMatrix[px][[17]]+DiagonalMatrix[px][[18]]}]

MatrixForm[leslie]

Eigenvalues[leslie]

es=Transpose[Eigensystem[leslie]][[1]]

The stable age distribution:

ev=es[[2]]/Sum[es[[2]][[i]],{i,1,18}]

Transpose[T2][[1]]/Sum[Transpose[T2][[1]][[i]],{i,1,18}]

Using the formulae for a Leslie matrix for the eigenvalue:

(1-Sum[lx[[i]] mx[[i]]/lam^i,{i,1,18}])//Factor

Solve[%==0,lam]

Using the formulae for a Leslie matrix for the right eigenvector:

Flatten[Join[Table[lx[[i]] (0.9771896634615088)^(-(i-1)),{i,1,17}],{lx[[18]](0.9771896634615088)^(-(18-2))(0.9771896634615088-(1-dx[[18]]))^(-1)}]]/(Sum[lx[[i]] (0.9771896634615088)^(-(i-1)),{i,1,17}]+lx[[18]](0.9771896634615088)^(-(18-2))(0.9771896634615088-(1-dx[[18]]))^(-1))

%-ev

The age class with the highest mi*ui:

Transpose[T4c][[1]]*ev

Projection using the leslie matrix:

T5=Transpose[{leslie.Transpose[T2][[1]],Transpose[T2][[2]]}]

The data from 1996:

Transpose[T2c][[1]]

Transpose[{Transpose[T5][[1]],Transpose[T2c][[1]]}]

BarChart[Transpose[{Transpose[T5][[1]],Transpose[T2c][[1]]}],AspectRatio->0.6,ChartLabels->{Transpose[T2c][[2]],None},
BarSpacing->{0,0.5},ChartStyle->{Gray,Black}]

From Statistics Canada (1996, female immigrants by age):

T6={{6956,"0"},{8210,"5"},{8303,"10"},{8569,"15"},
{11340,"20"},{14569,"25"},{14149,"30"},{11228,"35"},{7504,"40"},
{4541,"45"},{3157,"50"},{3081,"55"},{2568,"60"},
{1767,"65"},{1068,"70"},{552,"75"},{226,"80"},{110+22,"85+"}};

BarChart[Transpose[T6][[1]],AspectRatio->0.6,ChartLabels->Transpose[T6][[2]],BarSpacing->0.2,PlotRange->{Automatic,{-750,15000}}]

From Statistics Canada (1996, female emigrants by age):

T6a={{1059,"0"},{1723,"5"},{1830,"10"},{1538,"15"},
{1669,"20"},{3345,"25"},{3271,"30"},{2899,"35"},{2330,"40"},
{1384,"45"},{769,"50"},{455,"55"},{252,"60"},
{481,"65"},{300,"70"},{82,"75"},{45,"80"},{21+17,"85+"}};

BarChart[Transpose[T6a][[1]],AspectRatio->0.6,ChartLabels->Transpose[T6a][[2]],BarSpacing->0.2,PlotRange->{Automatic,{-200,4000}}]

From Statistics Canada (1996, returning female Canadians by age):

T6b={{718,"0"},{651,"5"},{596,"10"},{864,"15"},
{1748,"20"},{1729,"25"},{1311,"30"},{950,"35"},{727,"40"},
{476,"45"},{301,"50"},{314,"55"},{252,"60"},
{199,"65"},{130,"70"},{76,"75"},{47,"80"},{23+1,"85+"}};

BarChart[Transpose[T6b][[1]],AspectRatio->0.6,ChartLabels->Transpose[T6b][[2]],BarSpacing->0.2,PlotRange->{Automatic,{-200,4000}}]

Net migrants by age:

T6c=Transpose[{Transpose[T6][[1]]-Transpose[T6a][[1]]+Transpose[T6b][[1]],Transpose[T6][[2]]}]

BarChart[Transpose[T6c][[1]],AspectRatio->0.6,ChartLabels->Transpose[T6c][[2]],BarSpacing->0.2,PlotRange->{Automatic,{-750,15000}}]

Net female migrants:

Sum[Transpose[T6][[1]][[i]],{i,1,18}]-Sum[Transpose[T6a][[1]][[i]],{i,1,18}]+
Sum[Transpose[T6b][[1]][[i]],{i,1,18}]

Net migrants:

(208791-47230+22565)

Including net migration (note that we have guesstimated by multiplying
the number of migrants in 1996 by the length of the census):

BarChart[Transpose[{Transpose[T5][[1]]+5*Transpose[T6c][[1]],Transpose[T2c][[1]]}],AspectRatio->0.6,ChartLabels->{Transpose[T2c][[2]],None},
BarSpacing->{0,0.5},ChartStyle->{Gray,Black}]

(1) Projected female pop size, (2) Projected female pop size including migrants,
(3) 1991 female pop size, (4) 1996 female pop size

{Sum[Transpose[T5][[1]][[i]],{i,1,18}],
Sum[(Transpose[T5][[1]]+5*Transpose[T6c][[1]])[[i]],{i,1,18}],
Sum[Transpose[T2][[1]][[i]],{i,1,18}],
Sum[Transpose[T2c][[1]][[i]],{i,1,18}]}

BarChart[MatrixPower[leslie,20].Transpose[T2][[1]],AspectRatio->0.6,ChartLabels->Transpose[T2][[2]],BarSpacing->0.2,PlotRange->{Automatic,{-60000,1200000}}]

The reproductive values (left eigenvector) can be found by first transposing the Leslie matrix:

Transpose[Eigensystem[Transpose[leslie]]][[1]]

%[[2]]/%[[2,1]]

%.leslie-% 0.9771896634615126