FEBRUARY 1983 

TWO COMPUTER PROGRAMS TO PROJECT 

POPULATIONS WITH TIME-VARYING VITAL RATES 

Tim Gerrodette 

Daniel Goodman 

Jay Barlow 

NOAA-TM-NMFS-S WFC-28 

U.S. DEPARTMENT OF COMMERCE 

National Oceanic and Atmospheric Administration 

National Marine Fisheries Service 

Southwest Fisheries Center 
NOAA Technical Memorandum NMFS 

The National Oceanic and Atmospheric Administration (NOAA), organized in 1970, has evolved 

into an agency which establishes national policies an manages and conserves our oceanic, coastal, and 

atmospheric resources. An organizational element within NOAA, the Office of Fisheries is responsible for fisheries policy and the direction of the National Marine Fisheries Service (NMFS). 

In addition to its formal publications, the NMFS uses the NOAA Technical Memorandum series to 

issue informal scientific and technical publications when complete formal review and editorial processing are not appropriate or feasible. Documents within this series, however, reflect sound professional 

work and may be referenced in the formal scientific and technical literature. 
NOAA Technical Memorandum NMFS 

Thls TM series is used tor documentallon and tlmely communication of preliminary results, lmterim reports, or 

special purpose intormallon and have not recelved complete lormal review, editorial control, or detalled editing 

FEBRUARY 1983 

TWO COMPUTER PROGRAMS TO PROJECT 

POPULATIONS WITH TIME-VARYING VITAL RATES 

Tim Gerrodette and Jay Barlow 

Southwest Fisheries Center 

National Marine Fisheries Service 

La Jolla, CA 92038 

and 

Daniel Goodman 

Scripps Institution of Oceanography 

University of California, San Diego 

La Jolla, CA 92093 

N 0 A A-TM-N M FS-S W FC-2 8 

U.S. DEPARTMENT OF COMMERCE 

Malcolm Baldrige. Secretary 

National Oceanic and Atmospheric Administration 

John V. Byrne. Administrator 

National Marine Fisheries Service 

William G. Gordon, Assistant Administrator for Fisheries 
CONTENTS 

Page 

ABSTRACT ........................................................... 1 

INTRODUCTION ....................................................... 1 

PROGRAM DESCRIPTIONS ............................................... 1 

NUMERICAL EXAMPLE .................................................. 2 

LITERATURE CITED ................................................... 4 

APPENDIX 1 . User's guide to program SPP ........................... 14 

APPENDIX 2 . User's guide to program SLT ........................... 16 

APPENDIX 3 . Listing of program SPP ................................ 19 

APPENDIX 4 . Listing of program SLT ................................ 33 

i 
LIST OF TABLES 

Tab1 e Page 

1 Mean vital rates and covariance matrix of vital rates 

for the female population of the United States in 1940 

(from Sykes 1969). M refers to effective fecundity, P 

to survival rate, and numbers to 15-year age-classes.. ... 5 

2 Output of program SPP for the projection of the 1940 

U. S. female population for 7 time steps, using the 

vital rates given in Table 1 ............................. 6 

3 Output of program SLT for the projection of the 1940 

U. S. female population for 7 time steps. Sample size 

for the simulation was 10,000 trials. The results in 

section 2 of this table should be compared with the 

"predicted" values for time step 7 in Table 2 ............ 10 

ii 
TWO COMPUTER PROGRAMS TO PROJECT POPULATIONS 

WITH TIME-VARYING VITAL RATES 

Tim Gerrodette*, Daniel Goodman 

and Jay Barlow* 

Scripps Institution of Oceanography 

University of California, San Diego 

La Jolla, California 92093 

ABSTRACT 

Two computer programs which project age-structured 

populations using time-varying survival and fecundity rates are 

described and illustrated with a numerical example from human 

demography. Detailed instructions for the use of both programs 

are provided, and FORTRAN source code listings are given in 

appendices. 

INTRODUCTION 

This memorandum describes and illustrates the application of two computer 

programs developed to estimate future population size and age structure when 

survi Val and fecundi ty rates are variable. Due to our 1 ess-than-perfect abi 1 i ty to forecast future envi ronmental conditions , there wi 11 be some 

uncertainty about future population size, even if we have estimated the vital 

rates accurately under present conditions. We compute the future population 

vectors using a projection matrix (Leslie 1945). However, the elements of the 

Leslie matrix are assumed to be random variables. With variable vital rates, 

the population vector and total population size in the future are no longer 

fixed quantities, but variates with some (often unknown) distributions. The 

programs described in this memorandum compute the means and variances for a 

variety of demographic parameters when vital rates vary randomly. In 

particul ar, conf i dence i nterval s on popul ati on size and the probabi 1 i ty that 

the population will increase are computed. We will briefly describe the two 

programs, provide instructions for their operation, and give a numerical 

example of their application. The work was supported by NOAA Contract 

80-ABC-00207. 

PROGRAM DESCRIPTIONS 

FORTRAN program SPP (Stochastic Population Projection) is based on a 

formula derived by Sykes (1969). This program uses the third of three models 

presented by Sykes; in this model the elements of the Leslie matrix (the 

effective fecundity rates and the survival rates) are random variables, each 

*Present address: National Marine Fisheries Service, NOAA, Southwest 

Fisheries Center, P. 0. Box 271, La Jolla, CA 92038 
2 

with a specified mean and variance, and with specified covariances between all 

elements of the Leslie matrix. There is no autocovariance over time. 

Equations 19 and 20 of Sykes (1969) were used to compute the mean final 

population vector and the variance-covariance matrix for the final population 

vector, respec ti vel y . 

From these quantities the mean and variance of the total population can 

be computed. The expected value for the total population is the sum of the 

expected values of each age class, while the variance of the total population 

is the sum of the entries in the variance-covariance matrix. The distribution 

of the final total population is not normal, but confidence limits can be 

computed from the realized factor of increase. Mathematical details of the 

procedure are given in Gerrodette et a1 . (1983) .1 

FORTRAN program SLT (Stochastic Life Tab1 e simulator) simulates the 

growth of an age-structured population under fluctuating environmental 

conditions. At each time period, a new set of fecundity and survival rates 

(the elements of the Leslie matrix) are chosen and used to project the 

population. Each fecundity and survival rate is a normally distributed random 

variable with a mean, a variance, and covariances with every other fecundity 

and survival rate specified by the user. Distributions other than normal 

could also be used. The projection, starting from the same initial population 

vector, may be replicated a given number of times. From these results, the 

mean, variance, and covariances of the population vector are computed, 

together with statistics on a variety of other demographic parameters, 

including the distributions of the total population and the realized factor of 

increase. 

The two programs have been written in similar, easy-to-use formats. Each 

requires only an input file in which the data and control parameters are 

given; after that, each program runs automatically. Guides to the use of 

these two programs are given in Appendices 1 and 2. Source code listings are 

provided in Appendices 3 and 4. 

NUMERICAL EXAMPLE 

We illustrate the application of programs SPP and SLT by a projection of 

the female population of the United States in 1940, an example introduced by 

Keyfitz (1964) and used by Sykes (1969). Each age-class (and time step) in 

this example represents a span of 15 years. The mean vital rates and the 

covariance matrix for the vital rates over time are given in Table 1; we 

project the initial population vector through 7 time steps, a total of 105 

years, as Sykes does. The results of program SPP presented in Table 2 differ 

sl ightly from Sykes' ( 1969) results. After careful ly checking our results and 

after programming Sykes' formula in two different ways with identical answers, 

we believe that the differences can be ascribed to round-off errors owing to 

lGerrodette, T., D. Goodman, and J. P. Barlow. 

when vital rates vary randomly. In prep. 

1983. Population projection 
3 

differences in computational hardware. In addition, independent programming 

of Sykes' formula by T. Polacheck of the University of Oregon has given 

results identical to ours (T. Polacheck, Univ. of Oregon, Eugene, OR. pers. 

corn. 1. 

In Table 2, the expected population vector approaches the stable age 

distribution with each time step, as follows from Sykes' Equation 19. The 

covariance matrix for the population vector gives, on the diagonal, the 

variances of each age-class and, above the diagonal, the covariances of each 

age-class with every other age-class. The 95% confidence interval for the 

total population and for the factor of increase realized relative to the 

initial population are shown for each time step. After 30 years, for example, 

we can be 95% confident that the size of the population will be between 34.733 

and 60.276 million females, with a mean of 46.858 million. Because the 

distribution of final population size is skewed in a stochastic projection 

(Gerrodette et al. 19821, the mean will not, in general , lie midway between 

the lower and upper confidence limits. The probability that the total 

population will have increased over its initial value is also shown for each 

time step. Due to transients in the population vector, this probability 

decreases at first and then increases in this example. 

The results of the Monte Carlo simulation of the same population 

projection (program SLT) are given in Table 3. By comparing the results in 

Table 3, Section 2, with those of time step 7 in Table 2, we see that the 

results of the simulation (SLT) and the analytic solution (SPP) agree closely. 

As a check on the confidence limits for the total population computed by 

program SPP, program SLT tallies the proportion of final populations which 

fall above and below the computed upper and lower confidence limits. We 

expect that 2.5% of the cases should fall above the upper limit and 2.5% below 

the lower limit if the 95% confidence interval has been correctly estimated. 

The last part of Section 2 of Table 3 shows that 2.25% of the cases fell above 

the upper limit and 2.57% below the lower limit on this particular run. Other 

runs gave similar results. Program SLT also computes the proportion of cases 

in which the final population was greater than the initial population, and 

this answer (0.8664, Table 3) is close to the probability computed 

analytically by program SPP assuming that the realized factor of increase is 

normally distributed (0.8697, Table 2). Thus, if vital rates are variable but 

follow the distributions given in Table 1, we can say that the probability 

that the population will be larger 105 years in the future is 0.87, while the 

probability that it will be smaller is 0.13. 

The third section of Table 3 shows that the distributions of the vital 

rates achieved on this run were close to the "target" rates specified in Table 

1 (and repeated in the first parts of Sections 1 and 2 of Table 3). Section 4 

of Table 3 shows the results of a deterministic projection using the achieved 

mean values of the vital rates. For reasonably large simulation runs, Section 

4 should be similar to Section 1, which carries out a deterministic projection 

using the specified input values of vital rates. 
4 

LITERATURE CITED 

Keyfitz, N. 1964. The population projection as a matrix operator. Demography 

1~56-73. 

Leslie, P. H. 1945. On the use of matrices in certain population 

mathematics. Biometrika 33:183-212. 

Sykes, Z. M. 1969. Some stochastic versions of the matrix model for 

population dynamics. ---- J. Amer. Stat. Assn. 64:lll-130. 
Table 1. Mean vital rates and covariance matrix of vital rates for the female 

population of the United States in 1940 (from Sykes 196q). M refers 

to effective fecundity, P to survival rate, and numbers to If;-year 

age-cl asses. 

Mean 

MI 

MI n.?m7 .nmi 

M? R.68154 

M? n.imn 

PI n.wm 

P? n.07203 

P3 0.00000 
6 

Table 7. Output of program SPP for the projection of the 1940 U.S. female 

population for 7 time steps, using the vital rates given in Table 1. 

* * * * STOCHASTIC POPULATION PROJECTION * * * * 

GIVEN THE FDLLOWING INITIAL POPULATIDN ANn VITAL RATES - 

AGE CLASS POPULATION SIZE FECUNDITY (M) SURVIVAL RATE (P) 

1 14459.0n .32167 .98m 

2 15764.OO .C;P154 .97m 

3 11 346. no .12110 o. moon 

------------ 

TOTAL = 410m.no 

THE VITAL RATES - 

2 3 

. on1 88 n. noonn 

.nnim n.nonnc, 

- . nnm4 n. Omnn . 00n47 O. nnmn 

.om76 0. rwmn 

o . nnnnn O. nnnnr, 

THE EXPECTED POPULATION VECTOR AFTER 1 T’ME INTERVALS IS - 

ACE CLASS POPULATION SIZE 

1 16478.05 

7 14258.03 

3 14837.07 

TOTAL = 45523.14 

------------ 

WITH THE COVARIANCE MATRIX - 

1 3 3 

1 .13855Etfl8 .34595E+06 .64273E+n6 

3 .r)O0Or)Etl)l .00000E+nl .17707Et06 

7 . ~OOOOE~OI . Fin62m-05 . iom~tn6 

THE Q!iy CONFIDENCE LIMITS ARE - 

FflR THE FINAL POPULATION: .37614EtO5 TO .53432E+O!i 

FOR THE FACTOR OF INCREASE: .4159 TO 1. 3oin 

THE PRnRARILITY THAT THE POPULATION IS INCREASING JS .P657 

THE EXPECTED POPULATION VECTOR AFTER 3 TIME INTERVALS IS - 

AGE CLASS POPULATION SIZE 

I 16798.55 

> 16199.71 

3 13859.33 
7 

THE 957, CONFIDENCE LIMITS ARE - 

FOR THE FINAL POPULATION: .34733~+05 TO .~inm~+n~i 

FOR THE FACTOR OF INCREASE: .Ql% TO 1.2115 

THE PRORARILITY THAT THE POPULATION IS INCREASING IS .Pln8 

THE EXPECTEP POPULATION VECTOR AFTER 3 TIME INTERVALS IS - 

ACE CLASS POPULATION SIZE 

2 16565.P9 

I 18127.70 

3 15746. Fin 

------------ 

TOTAL = 50434.39 

WITH THE COVARIAFCE MATRIX - 

1 3 3 

1 .37783E+W .?1136EtO7 .114olEtO8 

z mnnmtni .15828~+08 .F;o~E+o~ 

3 .nnnnoE+ni .monE+ni .imi7~+ng 

THE 957, CONFTDENCE LIMITS ARE - 

FOR THE FINAL POPULATION: .~VE+OS TO .7zm~+n5 

FOR THE FACTOR OF TNCREASE: .9718 TO 1.7098 

THE PRORABILITY THAT THE POPULATION IS INCREASING IS .815O 

THE EXPECTEP POPULATION VECTOR AFTER 4 TIME INTERVALS IS - 

AGE CLASS POPULATION SIZE 

1 iqnx.71 

2 i787n.m 

3 16101.77 

TOTAL = 53998 78 

------------ 

THE 9% CONFIDENCE LIMITS ARE - 

FOR THE FINAL POPULATION: .31573E+O5 TO .80725E+O5 

FOR THE FACTOR OF INCREASE: .(I364 TO 1.3841 
8 

THE PROBABILITY THAT THE POPULATION IS INCREASING IS .8298 

THE EXPECTED POPULATION VECTOR AFTER S TIME INTERVALS IS - 

AGE CLASS POPULATION SIZE 

I 20249.75 

2 18761.75 

3 17370.05 

TOTAL = 56382.45 

------------ 

WITH THE COVARIANCE MATRIX - 

1 2 3 

1 .48931E+n8 .34973E+l18 .25746E+O8 

7 .m-mnnE+oi .wwE+n8 . I~ZWE+OR 

3 .nmnnE+ni .nmnnE+ni .?ti886~+08 

THE 95q CONFIDENCE LIMITS ARE - 

FOR THE FINAL POPULATION: .30847E+n5 TO .91465E+Wi 

FOR THE FACTOR OF INCREASE: .94n4 TO 1.1737 

THE PROBABILITY THAT THE POPULATION IS INCREASIVG IS .8436 

THE EXPECTED POPULATION VECTOR AFTER 6 TIME INTERVALS IS - 

ACE CLASS POPULATION SIZE 

1 71404.74 

3 19968.78 

3 18236.98 

------------ 

TOTAL = wFin9.f;~ 

WITH THE COVARIANCE MATRIX - 

1 2 3 

1 .F2527E+08 .35957E+n8 .33753E+08 

7 .nOnOOE+Ol .47604E+08 .?4O8C)E+O8 

3 .OnonnE+ni .nnnnnE+ni .33nx~+n8 

THF Q57, CONFIDENCE LIMJTS ARE - 

FOR THE FINAL POPULATION: .306lRE+n5 TI, .10131E+n6 

FPR THE FACTOR OF IF!CREbSE: ."5?3 TO 1.1624 

THE PRORARTLITY THAT THE POPULATION IS INCREASING IS .8575 

THE EXPECTED PnPIJLATION VECTOR AFTER 7 TIME INTERVALS IS - 

AGE CLASS POP[ILATION SIZE 

1 7270?.71\ 

7 31 106.73 

3 19409.77 
TOTAL = 63214.77 

THE 95:‘ COMFIDEFJCE LIMITS ARE - 

FOR THE FJNAL POPULATION: .~O~XE+M TO .II~E+M 

FOR THE FACTnR OF 1NCREASE: .OW1 TO 1.1549 

THE PRPRARILITY THAT THE POPULATION IS INCREASING IS .8697 
10 

Table 3. Output of program SLT for the projection of the 1990 U.S. female 

population for 7 time steps. Sample size for the simulation was 

10,nOO trials. The results in section 2 of this table should be 

compared with the "predicted" values for time step 7 in Table 2. 

* * * * * STOCHASTIC LIFE TARLE SIMULATION * * * * * 

1. INPIJT 

THE FOLLOWING VITAL PARAMETERS HAVE BEEN SPECIFIED: 

AGE BIRTH SURVIVAL REPRODUCTIVE 

CLASS RATE RATE VALUE DER/DM DER/DP 

1 .37167 .qnm 1 mooo . fim4 .412?7 

2 .68154 .Q72O.1 .74817 .51289 . n5~3 

3 .ui in n. nnmn .11431 .47057 o.0onnn 

ASYMPTOTIC FACTOR OF INCREASE I. 0594 

NET REPLACEMENT RATE 1.1098 

GENERATION TIME 1. Pl47 

VAR. IN AGE OF REPRODUCTION .36O1 

USING THESE VITAL RATES, THE 1NITIAL POPIJLATJON 1s PROJECTED 

FOR 7 TIME STEPS AS FOLLOWS: 

AGE INITIAL FINAL 

CLASS POPULATION PO P 11 L AT I 0 N 

1 14459. no 227n7.78 

2 15264. no ziim. 72 

3 11 346.00 19409.77 

F INAL/l NITlAL POPULATION I. 5393 

REALIZED FACTOR OF JNCREASE 1.0636 

2. STOCHASTIC POPULATION PROJECT1 ON 

THE FOLLOWING COVARIAVCE MATRIX HAS PEEN SPECIFIED: 

(FIRST 3 COLUMNS ARE BIRTH RATES, SECOND 3 ARE SURVIVAL RATES) 
11 

n.nnnnn n.nnmO n.nnnnn .now9 .nor147 o.nonno 

n.nnnnn n.onono o.oonnn n.nnnoo .on076 n.nOnno 

o.nnnoo n.noonn n.nnnno nmom o.nooon o.nnocIn 

USING THE MEAN VITAL RATES WITH THIS COVARIANCE STRUCTURE, 

THE INITIAL POPULATION IS PROJECTED FOR 7 TIME STEPS innon TIMES. 

AGE INITIAL (FIEAN) FINAL 

CLASS POPULATION POPULATION 

1 24459. vn 77682.34 

3 11346. no 19283.76 

7 15x4. nn 21085. 40 

CnVARIANCE MATRIX FOR FINAL POPULATION VECTOR 

MEAN VARIANCE 

FINAL/INITIAL POPULATION 1.5353 .75836 

REALIZED FACTOR OF INCREASE 1. R56l .no251 

FREOUENCY DISTR IRIITION FOR FACTOR OF INCREASE 

REALIZED 

LOWER BOUND FINAL/INITIAL LOWER ROUND FACTOR OF 

OF INTERVAL POPULATION OF INTERVAL INCREASE 

1 

2 

3 

4 

5 

6 

7 

8 

9 

11 

12 

13 

14 

1s 

16 

17 

18 

in 

< .mn 

. won 

.750o 

1. nom 

1.25~1 

1. !inon 

1.7500 

2. nnnn 

2. mn 

7. 5noo 

7.7500 

3. nnnn 

3. mo 

3. moo 

3.75n0 

4. nono 

4. xon 

.35O1) 

0 

262 

1803 

2095 

1779 

1296 

791 

457 

224 

107 

63 

31 

11 

4 

2 

in 

in64 

n 

< 325n 

.8750 

.850O 

.8750 

. ?eon 

.925O 

.wm 

.07m 

1. noon 

1.0250 

i.mn 

1. inno 

1.1500 

1.2000 

i.725n 

1. "500 

1.135n 

1.1751) 

n 

n 

130 

0 

5 

41 

374 

786 

1329 

1782 

1657 

1098 

543 

199 

71 

13 

2 

1970 
12 

19 4.5000 1 1.2500 0 

20 4.7500 0 1.2750 0 

PROPORTION OF PROJECTIONS WITH FACTOR OF INCREASE > 1.0 IS .8664 

PROPORTION OF FINAL/INITIAL POPULATION RATIOS BELOW .7410 IS .0257 

PROPORTION OF FINAL/INITIAL POPULATION RATIOS ABOVE 2.740 IS .0225 

3. STATISTICAL PROFILES OF DEMOGRAPHIC PARAMETERS USED IN STOCHASTIC 

POPULATION PROJECTIONS 

THESE STATISTICS ARE COMPUTED ON THE BASIS OF 70,000 VALUES 

( 7 TIME STEPS, 10000 REPETITIONS) 

AGE BIRTH RATE SURVIVAL RATE REPRODUCTIVE VALUE DER/DM DER/DP 

CLASS MEAN VAR MEAN VAR MEAN VAR MEAN MEAN 

MEAN VARIANCE 

ASYMPTOTIC FACTOR OF INCREASE 1.0564 .nip17 

NET REPLACEMENT RATE 1.1086 .o6m 

GENERATION TIME 1.8250 .OO627 

VAR. IN AGE OF REPRODUCTION .3558 .00060 

***** WARNING: FACTOR OF INCREASE COULD NOT BE COMPUTED 14 TIMES 

THE FOLLOWING COVARIANCE STRUCTURE WAS ACHIEVED: 

4. DETERMINISTIC POPULATION PROJECTION 

THE ACHIEVED MEAN VITAL RATES PRODUCE THE FOLLOWING ASSOCIATED 

PARAMETERS: 

AGE BIRTH SURVI VAL REPRODUCTIVE 

CLASS RATE RATE VALUE DER/DM DER/DP 
13 

ASYMPTOTIC FACTOR OF INCREASE 1.0589 

NET REPLACEMENT RATE I. 1089 

GENERATION TIME 1.8156 

VAR. IN ACE OF REPRODIICTION .36n7 

USING THESE ACHIEVE0 MEAW VITAL RATES, THE INITIAL POPULATION 

IS PROJECTED FOR 7 TIME STEPS AS FOLLOWS: 

AGE INITIAL FINAL 

CLASS POPULATION POPULATION 

1 14459. no ?2m7.48 

3 11346. fV 19275.70 

2 15264.00 mi 3.83 

F 1 NAL/ I NIT I AL POPIJLATI ON 1.5325 

REALIZED FACTOR OF JRCREASE 1. ~29 
APPENDIX 1. User's guide to program SPP. 

A. General 

The program was written in Microsoft F80 FORTRAN. This FORTRAN is 

compatible with almost any machine except for machine specificity of the 

I/O (READ, WRITE, and OPEN statements). The random number generator 

(subroutine RANDM) may a1 so be rnachine-specific. 

R. Input 

The input file, called '1NPUT.SPP' , contains the control parameters 

and data in the following order and format: 

Record 1: N 

Record 2: NT 

Record 3: ICONF 

Record 4: IOPl 

Record 5: IOP? 

Record 6: IPP3 

Record 7: IOP4 

Record 8: EMX 

Record 9: PX 

Record 10: X 

Record 11: CllVAR 

The number of age classes in the population 

to be projected. As currently dimensioned, the 

program handles up to 75 age classes (I?) 

The number of time steps the population is to be 

projected (I?) 

The percentage confidence limits desired; (1.0- 

alpha error) x 100; thus if 95 is specified here, 

the 955 confidence limits will be computed (I?) 

An option (implemented by IllP1>0) which allows 

the initial population vector to be the stable 

age vector; otherwise (IOPl=ll) the initial 

population vector is specified in the input file 

- see Record 10 below (11: 

An option (implemented by IOP?>n) which allows a 

simp1 ified covariance structure to be specified; 

otherwise ( IOP2=0) covariances between each 

fecundity and survival rate must be specified 

in the input file - see Records 11-14 below (11) 

An output option (implemented by IOPR>O> which 

lists the input data (Records 8-11); otherwise 

(IOP3=ll) only the results are output (11) 

An output option (implemented by IOPfi>O) which 

displays the results of the computations after 

each time step; otherwise (IPP4=0) only the 

final results after NT time steps are shown (11) 

Effective fecundity rate vector; the first row of 

the Leslie matrix (NF10.6, where N is Record 1) 

Survival rate vector; the subdiagonal of the 

Leslie matrix (NF10.6, where N is Record 1) 

If IOPI=n, the initial population vector (NF10.2); 

if IOPl>O, the total initial population (F10.2) 

lf IllP?=O, the elements of the covariance matrix 

of the fecundity and survival rates. COVAR is 

visualized as an upper triangular matrix, but 

entered as a vector by reading from left to right 

across each row, thus: 
15 

M1 M2 P1 

M1 123 

M2 56 

P1 8 

P2 

If I0P2>OY C0VAR conta 

fecundity and survi Val 

in the example above) 

P2 

4 

7 

9 

10 (N( 2N+l)Fln.6) 

ns the variances of the 

rates (elements 1,5,8,10 

2NF10.6) 

(Records 12-14 are read only if IOP2>0) 

Record 12: C1 The correlation among all fecundity rates (Flf1.6) 

Record 13: C2 The correlation between a1 1 fecundity and a1 1 

Record 14: C3 The correlation among a1 1 survival rates (F10.6) 

survival rates (Fln.6) 

The following diagram will clarify the use of IOP2, 

the simplified covariance option: 

M1 M2 . . . Mn P1 P2 . . . Pn 

M1 V1 

M2 v2 c1 . c2 

Mn 

P1 

P2 

Pn 

Vn ................ 

Vn+2 C3 

Vn+l 

V2n 

where the V's are the variances specified in 

Record 11. 

C. Execution 

Once the input file is specified, the program is executed by 

the command SPP. The input file and program SPP must be on the 

'current disk drive.' 

D. output 

The output is written to the console screen and to a diskfile 

called '0IlTPUT.SPP'. 

E. Other comments 

Caution should be exercised in specifying the variances and 

covariances of the vital rates. Large variances may cause the 

vital rates to assume unreasonable values (e.g., a survival rate 

greater than 1.0 or a negative fecundity rate). The program does 

not detect thi s. Variances must be positive , whi 1 e covariances 

may be either negative or positive. Zeros are permissible. 
16 

APPENDIX 2. User's guide to program SLT. 

A. General 

The program was written in Microsoft F80 FORTRAN. This FORTRAN is 

compatible with almost any machine except for machine specificity of the 

I/O (READ, WRITE, and OPEN statements). The random number generator 

(subroutine RANDM) may a1 so be machine-specific. 

6. Input 

The input file, called '1NPUT.SLT' , contains the control parameters 

and data in the following order and format: 

Record 1: 

Record 2: 

Record 3: 

Record 4: 

Record 5: 

Record 6: 

Record 7: 

Record 8: 

Record 9: 

Record 10: 

Record 11: 

Record 12: 

Record 13: 

Record 14: 

N 

NT 

NR 

IOPl 

IOP2 

LU N 

ISEED 

KK 

WIDTH1 

BLI. 

WIDTH2 

RL2 

RL 

R IJ 

The number of age classes in the population 

to be projected. As currently dimensioned, the 

program handles up to 20 age classes (12) The number of time steps the population is to be 

projected (I?) 

The number of times the simulation is to be 

repeated; the "sample size" of the stochastic 

simulation (15) 

An option (implemented by IOP1>0) which allows 

the initial population vector to be the stable 

age vector; otherwise (IOP1=0) the initial 

population vector is specified in the input file 

- see Record 17 below (11) 

An option (implemented by IOP2>0) which allows a 

simplified covariance structure to be specified; 

otherwise ( IOP2=0) covariances between each 

fecundity and survival rate must be specified 

in the input file - see Records 18-21 below (11) 

The logical unit number to which the output is 

directed during execution (e.g., console), other 

than 1 or 6, the input and output LUNs (11) 

The seed for the random number generator; a 

positive integer in the range 1-32767 (IS) The number of intervals in the frequency 

distribution to be computed (12) 

The width of the interval for the frequency 

distribution of the factor of increase (F10.6) 

The lower bound of the frequency distribution 

of the factor of increase (F10.6) 

The width of the interval for the frequency 

distribution of the final/initial population 

ratio (F10.6) I 

The lower bound of the frequency distribution 

of the final/initial population ratio (F10.6) 

Estimated lower final/initial population ratio 

(Flfl.6) 

Estimated upper final/initial population ratio 

(F10.6) 
17 

Record 15: RX Effective fecundity rate vector; the first row of 

Record 16: PX Survival rate vector; the subdiagonal of the 

Record 17: RNX If IOPl=O, the initial population vector (NFln.7); 

Record 18: S If 10P?=O, the elements of the covariance matrix 

the Leslie matrix (FJF10.6, where N is Record 1) 

Leslie matrix (NFlO.6, where N is Record 1) 

if IOPl>O, the total initial population (Fln.2) 

of the fecundity and survival rates. S is visuali zed as an upper triangul ar matrix, but entered 

as a vector by reading from left to right across 

each row, thus: 

M1 M2 P1 P2 

M1 1234 

M? 567 

P1 89 

P2 10 (N( 2N+1 )F10.6) 

If IT)P2>T), S contains the variances of the 

fecundity and survival rates (elements l,E;,R,lO 

in the example above) (7NF10.6) 

(Records 19-21 are read only if 10P2>0) 

Record 19: C1 The correlation among all fecundity rates (F10.6) 

Record 20: C? The correlation between all fecundity and all 

Record 21: C3 The correlation among all survival rates (Fln.6) 

survival rates (FlO.6) 

The following diagram will clarify the use of IOP2, 

the simp1 ified covariance option: 

M1 M2.. . F.ln P1 P2. . Pn 

M1 V1 

M? v2 c1 . c2 

hl 

P1 

P2 

Pn 

. Vn ................ Vnt2 C3 

Vn+l 

. 

V 2n 

where the V's are the variances specified in 

Record 18. 

C. Execution 

Once the input file is specified, the program is executed by 

the command SLT. The input file and program SLT must be on the 

'current disk drive.' 
18 

n. Output 

The output is written to the device specified by LUN and to a 

diskfile called 'OUTPUT.SLT'. The output is in 4 parts. First, the 

vital rates specified as input are listed, together with associated 

demographic parameters such as reproductive value and asymptotic 

factor of increase. The population structure which results from 

projecting with these vital rates is shown. 

are shown. The results are summarized in terms of the mean and 

variance for each parameter. The frequency distributions of the 

ratio of the final to the initial population and of the realized 

factor of increased are given. Also shown are the proportions 

of the projections which fall below and above the estimated 

confidence limits given on input, and the proportion of the 

final population ratios which are larger than the initial popu1 at i on. 

parameters achieved in the projections is displayed. These 

statistics indicate how closely the simulation was able to match 

the specified vital rates and their covariance structure. If all 

has gone well, the values in this section should closely match the 

values in Section 1 of the output. Because the vital rates are 

randomly generated, it is possible that the algorithm to compute 

the asymptotic factor of increase will not converge. A message 

will appear if this has happened. 

Fourth, a deterministic population projection is carried out 

using the mean values actually obtained in the stochastic projection 

These results allow a direct comparison with the results in Section 2. 

Second, the results of the stochastic population projections 

Third, a statistical summary of the vital rates and associated 

E. Other comments 

In specifying the vital rates, variances must be positive, 

while covariances may be either negative or positive. Zeros are 

permissible. [In1 ess these variances and covariances are careful ly 

specified, the vital rates may assume unrealistic values such as 

negative survival rates. The program truncates the realized 

distribution so that the survival rates lie between 0 and 1 

(inclusive) and the fecundities are equal to or greater than 0. 

This truncation may cause the achieved means to deviate systematically from the specified values. 
19 

APPENDIX 3. Listing of program SPP. 

PROGRAM SPP - - STOCHASTIC POPULATION PROJECTION 

C 

C...THIS PROGRAM ESTIMATES THE EXPECTEP POPULATIflN VECTOR ANP ASSOCIATED 

c STATISTICS AFTER ITERATIVE APPLICATIONS OF A STOCHASTIC LESLIE MATRIX, 

C BASEII 0N THE SYKES ALGORITHM (1969,MODEL 3). THE CONFIDENCE LIMITS FOR 

C THE FACTOR OF INCREASE AND THE PRORARILITY THAT THE POPULATION IS 

C INCREASING ARE COMPUTED RY A FORMULA DEVELOPED FOR THIS PROGRAM. 

C.. .INPUT FILE: 

C RECORP 1: N=NUMBER OF AGE CLASSES (I?) 

C RECORD 7: NT=NUMBER OF TIME STEPS TO RE PROJECTED (I?) 

C RECORD 3: ICONF=CQNFIDENCE LIMITS @ESIREP (12) 

c RECORD 4: IOPI=OPTION FOR INITIAL STABLE ACE DISTRIBUTJON (11) 

C RECORD 5: IOP2=OPTION FOR SIMPLIFIED COVARIANCE STRUCTURE (11) 

c RECORD 6: IOP3=OPTION FOR LISTING OF INPUT DATA (11) 

c RECORP 7: IOP4=OPTION FOR OUTPUT AFTER EACH TIME STEP (11) 

c RECORD 8: EMX=VECTOR OF FECUNDITIES ( NF 1 0.6) 

C RECORD 9: PX=VECTOR OF SURVIVAL RATES (NF 10.6) 

c RECORP 10: X=INITIAL POPULATION VECTOR (NFlQ.7) 

C RECORD 11: COVAR=COVARIANCE MATRIX FOR VITAL RATES. COVAR IS VISUALIZEP 

C AS AN IJPPER TRIANGULAR MATRIX, RIIT ENTERED AS A VECTOR RY 

c READING FROM LEFT TO RIGHT ACROSS EACU ROW, E.G., 

C M1 M2 PI P? 

C MI 1 2 3 4, 

c M7 567 

C P1 89 (N (2N+l)F10.6) 

C P2 IO 

c IF THE SIMPLlFIED COVARIANCE OPTION IS USEP, COVAR IS 

C ENTEREP IN A DIFFERENT FASHION. REFER TO FILE SPPGUIDE. 

C...OUTPUT IS WRITTEN TO THE CRT (LUN=l) APJP TO DISKFILE 'OUTPUT.SPP'. 

C.. . REOUIRED SURROUTINES : SYKES ,LOAP, PTOL, MMNT, ERNWT, STABL, LMPHI , PRJCT,PS? , C PU,CLPOP,MOMENT,PMPY,PSlJRyPILPyP~SDy~R~NVyClJMPR. 

C THE VAR IANCE-COVAR IANCE STRUCTURE OF THE EXPECTEP POPULATION VECTOR IS 

DIMENSION X(75) ,EMX(25) ,PX(25) ,V(25,75) ,COVAR(1775) 

DIMENSION c~(inn),w(ioo),~(~7s) 

C 

C. .STEP ONE - - INPUT PARAMETERS AND COMPUTE STARLF: AGE PISTRIRIJTION 

CALL OPEN (10,'INPUT SPP' ,O) 

CALL OPEN (2,'nUTPUT SPP' ,n) 

CALL LOAII (N,NT,ICONF,IOP1,IOP7,1@P3,IOP~,EMX,PX,XyCOVAR) 

N7=N*2 

IF (10~1) 2~i,25,in 

10 POP=X(l) 

CALL PTOL (PX,X,N) 

CALL LMPHI (X,EMX,W,N) 

CALL MMNT (W,RZyGT,S2,N) 

CALL ERNWT (W ,RZ ,GT, S3, ER ,N) 

CALL STABL (ER , X , X , !I) 

DO 15 I=l,N 

15 X(I)=X(I)*POP 

C 

C...STEP TWO - - INITIALIZE PARAMETERS 
20 

25 ENO=0.O 

DO 311 I=l,W 

ENO=ENO+X(I) 

n0 30 J=l,N 

30 V( I , J)=0. 0 

c 

C...STEP THREE - - WRITE HEAPING AND ECHO INPUT PATA 

IF (IOP3) 70,70,35 

no An 1=1,N 

WRITE (2,705) EN0 

IF (PX(N)) 45,45,42 

35 WRITE (2,700) 

40 WR I TE (2 , 701 I , x ( I , EMX ( I , PX ( I 

47 WRITE (2 , 707) 

45 WRITE (2,303) (J,J=~,N),(J,J=~,N) 

JS T P=0 

P0 60 I=1,N7 

J STRT=J STPt 1 

(1 S TP=J STRT+N 7- I 

w( I )=n.p 

K=I-1 

IF (K) 50,50,55 

GO TO 60 

50 WRITE (3,713) I ,(C(~VAR(J) ,J=JSTRT,JSTP) 

55 IF (I-N) 56,56,57 

56 WRITE (2,713) I ,(W(J) ,~1~1,K),(~~VAR(~1),J~~1STRT,~1STP) 

57 IP=I-N 

GO TO Fin 

WRITE (3,714) IN, (W( J) ,J=1 ,K) , (COVAR(J) ,J=JSTRT,JSTP) 

60 CPNTI NUE 

c 

C...STEP FOUR - - PROJECT COVARIANCES IISING THE SYKES ALGORITHM 

70 DO 180 IT=l,NT 

CALL SYKES (F1 , X , V,35 , EMX ,PX,COVAR ,C ,W) 

C 

C...STEP FIVE - - PROJECT THE MEAN POPULATION VECTOR 

CALL PRJCT (X,PX,EMX,&) 

IF (IOP~ i45,145,15n 

145 IF (IT-NT) itwyiw,i5o 

c 

C...STEP SIX - - CPMPIJTE CONFIDENCE LIMITS 

c 

C.. .STEP SEVEN - - 0(ITPUT RESULTS 

150 CALL CLPOP (X,V,7~,N,1T,1~ONF,~N(1,E~T,W,C,CC,CLL,CL~l,P~PL,P~P~1,PR) 

WRITE (1,307) IT 

WRITE (2,3n7) IT 

PO 165 i=i,r! 

WRlTE (1,7nl) I,X(I) 

165 WRITE (2,701) I ,X( I) 

WRJTE (1,3(75) ENT 

WRITE (?,?nFi) ENT 

WRITE (1,710) (J,J=l,N) 

DO 171) I=l,V 

WR I TE ( 3 , 71 0) (J , J=1 , PI) 
21 

W( 1)=0.0 

K=I-1 

IF-(K) 168,168,169 

1m WRITE (1,311) I,(V(I,J),J=I,N) 

WRITE (2,211) I ,(V( I ,J) ,J=1 ,v) GO TO 170 

WRITE (7,211) I, (W( J) ,J=1 ,K) ,( V( I , J) ,J=I ,N) 

I70 CONTINUE 

IF (PR) 174,175,175 

I74 WRITE (1,708) 

WRITE (2,708) 

GO TO lP0 

175 WRITE (1,715) ICONF,P0PLYPOPU,CLL,CLU,PR 

WRITE (2,715) ICONF,POPL,POPU,CLL,CLU,PR 

180 CONTIN\IE 

169 WRITE (1,711) I,(W(J) ,J=l,K) ,(V( 1,J) ,J=I,N) 

c 

C. . . OUTPUT FPRMATS 

700 FORMAT (lOX,'* * * * STOCHASTIC POPULATION PROJECTION * * * *I// 

#lX,'GIVEN THE FOLLOWING INITIAL POPULATION AND VITAL RATES -I/ 

@X,'ACE CLASS',3XY'POPULATION SIZE' ,3X,'FEC!lNDITY (11)' ,3X,'SURVIVA 

#L RATE (P)') 

701 FORMAT (19,F19.7,7F16.5) 

707 FORMAT (/SX,'(YOU HAVE SPECIFIED A NON-ZERO SURVIVAL RATE FOR', 

if' THE LAST ACE CLASS.'/fiX,'IF IT IS TRULY THE LAST AGE CLASS, NO', 

8' SIIRVIVnRS ARE PERMJTTED.) I) 

FORMAT (/lX,'WITH THE FOLLOWING COVARIANCE STRUCTURE AMONG THE' , f' VITAL RATES -'/5(4X,lOIlO/)) 

703 

205 FORMAT (16XY12('-')/6X,'TOTAL =',Fl5.2) 

207 FORMAT (//lX,'THE EXPECTEP POPULATION VECTOR AFTER' ,13,' TIME', 

208 

#' INTERVALS IS -'/4X,'AGE CLASS',3X9'POPULATION SIZE') 

FORMAT (/IX,'YOUR SPECIFIED COVARIAKE STRIJCTURE AMONG THE VITAL' , f' RATES HAS LEP Tn '/lX,'AN AFIAMOLOIJS RESULT FOR THE EXPECTEP', 

ff' VARIANCE IN THE FACTOR OF INCREASE.'/lX,'FURTHER STATISTICS' , #' CANNOT BE COMPUTED. '/) 

210 FORMAT (/IX,'WITH THE COVARIANCE MATRIX -i/qix,inIii/)) 

21 1 

213 FORMAT (I M(',1?,')',5(lOFlO.!i/6X)) 

214 FORMAT (' P(',17,')',5(10FlO.5/6X)) 

215 FORMAT (/lX,'THE 'I?'% CONFIDENCE LIMITS ARE -I 

FCIRMAT ( 14,l X, 3( lnEl1.5/5X) ) 

#/5X,'FOR THE FINAL POPULATION:',E14.5,' TO',E12.5 

B/SX,'FOR THE FACTOR OF INCREASE:' ,F12.4,' TO' ,F12.4 

#//lX,'THE PROBABILITY THAT THE POPULATION IS INCREASING IS',F7.4/) 

CALL EXIT 

END 

C*****SUBROUTINES FOR PROGRAM SPP 

C 

C.....RETURNS COVARIANCE MATRIX V OF AN N-ORPER POPULATION VECTOR X 

SUBROUTINE SYKES (N, X ,V, MR , EMX, PX ,COVAR ,C ,W) 

C WHOSE DYNAMICS ARE GOVERNED BY THE FECUNDITY VECTOR EMX ANI, THE 

C SllRVIVAL RATE VECTOR PX WITH COVARIANCES AMONG VITAL RATES 

C SPECIFIED IN COVAR. NR IS THE ROW DIMENSION OF V IN THE CALLING 

c PROGRAM; W AND C ARE WORK ARRAYS WITH DIMENSIONS AT LEAST N AND 
32 

C.....N(N+l)/Z, RESPECTIVELY. THE PIMENSION OF COVAR IS N(2N+1). 

c ..... REFERENCE: Z.M.SYKES (iw), J.AM.STAT.ASSN. miii-i3n, EO. 20. 

DIMENSION V(NR,1>,X(l),EMX(l),PX(1),C~VAR(l),C(l),W(l) 

SCRPT( I , J,N)=( I-1)*N- ( (1-3)*1)/2-I+J 

N2=W*2 

suMi =n. n 

DQ 9n I=I,N 

 SUM^=^. n 

DO 8n J=I,N 

IF (LJ) Fio,fin,Fi5 

50 I J=SCRPT( I, 3, N2) 

GO TO 60 

55 I J=SCRPT(J , 1, N2) 

60 SUM~=SUM~+C~VAR(IJ)*(V(I,J)+X(I)*X(J)) 

IF (1-1) 80,80,65 

SUM 2= S UM 2+C WAR ( I J ) * ( V ( I - 1 , J ) +X ( I - 1 ) *X ( J ) ) 

65 IJ=SCRPT(J,I+N-l,N2) 

IF (~-1) m,~n,70 

70 IF (I-J) 75,75,nn 

75 I J=SCRPT( I-l+N , 3-l+N ,N2) 

K=SCRPT( I ,J,N) 

C( K)=COVAR (I J)* (V ( 1-1, J-1)+X ( I -l)*X (J-1) ) 

C( I )=SUM2 

C( l)=SIJMl 

80 CONTINUE 

90 CONTINUE 

no 110 J=I,N 

nn ino I=~,N 

inn W( I)=v( I,J) 

110 V( I ,J)=w( I) 

CALL PRJCT (W,PX,EMX,rJ) 

DO 110 I=l,N 

DO 130 J=l,N 

no im I=~,N 

DO 130 I=I,FJ 

no 140 1=i,~1 

IJ=SCRPT( I,J,N) 

inn V(J, I >=v( I, J) 

130 W( l)=V(J,I) 

CALL PRJCT (W,PX,EMX,N) 

130 V(J,I)=W(l) 

PO 140 J=I,N 

V( I ,J)=V( I ,J)+C( IJ) 

RETURN 

END 

SIJRROUTINE LOAD (N,NT,1CONF,1OP1,1OP2,1OP3,1OP4,EMX,PX,X,C0VAR) 

C.. . . .REAPS IWPUT DATA FOR PROGRAM SPP, INCLUDING REPACKING OF COVARIANCE 

C.. . . .MATRIX GIVEN SIMPLIFIED COVARIANCE STRUCTIJRE. 

DIMENSION X(l) ,EMX(l) ,PX(l) ,COVAR(l) 

1 FORMAT (11) 

2 FORMAT (I?) 

3 FORMAT (8F10.2) 

4 FORMAT (4F10.6) 
23 

READ (10,2) N 

N2=N*2 

READ (10,2) NT 

READ (lo,?) ICONF 

READ (10,l) IOPl 

READ (10,l) IOP3 

READ (10,l) IOPn 

READ (10,4) (EMX(I),I=l,N) 

REAP (lO,n) (PX( I) , J=l ,N) 

IF (IOP1) 15,15,10 

10 READ (10,3) X(l) 

GO TO 20 

15 READ (10,3) (X(I),I=l,N) 

70 IF (IOP2) 25,75,30 

75 READ (10,4) (COVAR(I),I=l,NCV) 

30 READ (l0,d) (COVAR(I),I=l,W2) 

NCV=N*(N?+I) 

READ (10~1) 10~2 

C.. . . .COVARIANCE VECTOR 

GO TO 2a-m 

READ (lO,4) C1 

RFAD (10,4) Cil 

READ (10,4) C3 

IF (ARS(Cl)-l.) 40,40,9fJ 

50 IF (ARS(C3)-1.) 6OY6c1,90 

40 IF (ABS(C2)-1.) 50,50,90 

60 DO 70 I=l,N2 

70 CONTINUE 

8r\ WRITE (1,RS) 

IF (COVAR(I)) 80,7n,70 

61) TO inno 

89 FORMAT (/I ERROR -- VARIANCES IN COVARIANCE MATRIX CANNOT'/ 

#9X'HAVE NEGATIVE VALUES I) 

GO TO 3200 

90 WRITE (1,99) 

99 FORMAT (/I ERROR -- ELEMENTS OF CORRELATION MATRIX CARNOT'/ 

P9X'HAVE ABSOLUTE VALUE GREATER THAN 1.0') 

GO TO 2200 

C.....REPACK 

1000 MIn=NCV-N*(N+l)/il 

NDX=NCV+l 

K=NCV 

I=N2 

M=l 

DO 1200 NNN=l,NCV 

NDX=NDX-1 

IF (NDX-K) 1100,l 160, 1100 

1100 IF (NDX-MID) iiio,iiio,iizn 

1110 IF (NDX-L) 1150,1140,1140 

1130 COVAR( NDX)=C3*SQRT(COVAR( I )*COVAR(J) ) 

1 140 C OVAR ( NDX) =C 2*SORT ( COVAR ( I ) "COVAR (3 ) ) 

1 1 50 C OVAR ( ND X ) =C 1 *SORT ( C OVA R ( I ) *COVAR ( J ) ) 

GO TO izno 

GO TI, 1200 
24 

GO TO izon 

iim COVAR(NDX)=COVAR(I) 

L=K-N 

M=M+1 

K=K-M 

I=I-1 

J=N?+l 

1200 J=J-1 

C.....ZERO-MEAN CASES 

3000 DO 3100 I=~,N 

mi0 IF (px(I)-.noooooi) mo,znzn,~m 

znm IR=I+N 

IF (EMX(I )-.OOOOOOI) z030,~o3n,znin 

GO TO 3040 

IF (J-IR) ~05n,3n5n,20m 

7030 IR=I 

3040 DO 7070 J=l ,N? 

2050 IJ=(J-l )*N?-( (J-3)*J)/?-J+IR 

7060 IJ=( IR-l)*N2-( (IR-3)*IR)/?-IR+J 

707n COVAR ( I J)=0. O 

3100 CONTINUE 

2300 RETURN 

ENP 

GO TO 7070 

IF (IR-LE) znin,mn,~in~ 

SIIRR(IUTTF\’E PRJCT( ENX,PX, EMX ,F) 

C*****PPVAVCES A POPIlLATInN VECTOR ONE TIME IlNIT FORWARD RY APPLICATION 

C OF A LESLIE MATRIX DEFINED RY N-ORPER VECTORS OF FECUNDITY AND 

C.. . . .SIlRVIVAL RATES. 

DIMENSION EFJX(1) ,PX(l),EMX(l) 

SIIM=0. 0 

I X=N+I 

I XMl=N 

DO 1050 T=Z,N 

IX=I x-1 

IXMI=IXM1-1 

SIlM=SllMtENX( I X)*EMX( IX) 

ENX ( I X) =ENX ( I XM1) *PX ( I XM1) 

1050 CONTINUE 

ENX ( ~)=S~IM+ENX (1 )*EMX (1 ) 

RETURN 

END 

SIIRROUTINE PTOL(PX,FLX,N) 

DIMENSION PX(1) ,ELX(l) 

ELX( 1 )=1. 0 

ELAST=l.O 

IM1=0 

no l05O I=?,FI 

TMl=IMl+I 

ELAST=ELAST*PX(IMI) 

EI.X( I )=FLAST 

C*****CPNVERTS AN N-ORDER VECTOR OF SURVIVAL RATES TO SIIRVIVORSHIPS. 
25 

1O5n CONTINUE 

RETllRN 

EN? 

SUBROUTINE LMPHI (ELX, EMX, PHI X,N) 

C*****CONVERTS AN N-ORDER VECTOR OF SURVIVORSHIPS AND A VECTOR OF 

C.. . . .FECUNDITIES TO A NET MATERNITY FUNCTION. 

DIMENSION ELX(1) ,EMX(l) ,PHIX(l) 

DO in5n I=I,FI 

PHI X( I )=ELX( I )*EMX ( I ) 

in50 CONTINUE 

RETURN 

END 

SlIRROUTINE MMF!T( PHIX,RZ ,T, S2,N) 

C*****CALCULATES THE NET REPLACEMENT RATE, THE MEAN T, AND THE VARIANCE S? 

C IN AGES OF REPRODUCTION FOR A COHORT. 

DIMENSION PHIX(1) 

RZ=fl. (I 

XlPH I=O. 0 

X?PH I=n. n 

DO ion0 I=~,N 

PHI=PH I X( I ) 

FI=FLOAT( I) 

RZ=RZtPH I 

XlPHI=XlPHItPHI*FI 

X?PHI=X2PHItPHI*FI*FI 

T=XlPHI/RZ 

RETURN 

END 

lPnO CONTINUE 

S?=XZPHI/RZ-T*T 

SURROUTINE ERNWT (PHI X ,RZ ,T, S?, ER ,N) 

C*****CALCIILATES THE ASYMPTOTIC FACTOR OF INCREASE (E RAISED TO THE POWER 

C OF THE MALTHUSIAW PARAMETER) RY NEWTON'S METHOD. THE INITIAL VALUE 

C.....IS CALCULATED FROM INPUT MOMENTS OF THE NET MATERNITY FUNCTION. 

DIMENSION PHIX(1) 

ER=EXP ( (T/S?)*(RZ** (S?/T**2)-1 .O) ) 

CT=T). OO01 

ERX=l. 0 

DFDR=n.O 

ERX=ERX*ER 

TRM=PHIX( I)/ERX 

F R=FR+TRM 

PFPR=DFDR -F LOAT ( I )*TRM 

DELT=FR/DFPR 

ABFR=ARS(FR) 

DO ~nnn J=I,?O 

FR=-1.0 

DO inw I=I,N 

ln5O CONTINUE 

IF(ARFR-CT) inm,inm,inw inm ARDLT=-ARFR/DFPR 
26 

IF (ABDLT-CT) ~05o,~n50,io9n 

1090 ER=ER-DELT 

70110 CONTINUE 

9nO0 FORMbT (/5Xy I ***** WARNING: NO CONVERGENCE IN SURROUTINE ERNWT'/) 

WRITE (2,9000) 

ER=-100.0 

RETURN 

END 

2050 ER=ER-DELT 

SUBROUTINE STABL (ER, ELX,CX ,N) 

C*****CALCULATES THE STABLE AGE DISTRIBUTION CX , GIVEN THE ASYMPTOTIC 

C*****FACTOR OF INCREASE AND THE SURVIVORSHIP VECTPR OF ORDER N. 

OIMENSION ELX(1) ,CX(l) 

SUM=O. 0 

EMRX=l .0 

EMR=l.O/ER 

Do 1000 I=l,N 

EMRX=EMRX*EMR 

CX( I)=ELX( I)*EMRX 

1000 SIJM=SUMtCX ( I ) 

DO 1100 I=l,N 

1100 CX( I)=CX( I)/SlIM 

RETURN 

END 

SIJRROIJT I NE MOMENT (N , C , NTERMS) 

C.....COMPUTES THE VECTOR OF COEFFICIENTS C OF THE N-TH MOMENT OF A 

C.. . . .NORMALLY DISTRIRUTED VARIABLE (THE N-TH DERIVATIVE OF THE MOMENT 

C.....GENERATING FUNCTION). 

C.....THE N9RMAL VARIABLE, VECTOR C WILL CONTAIN COEFFICIENTS OF A 

C.....EXPONENTS 0,1,?,3,... AND WHERE NTERMS IS THE NUMBER OF TERMS 

C.....IN THE POLYNOMIAL (DIMENSION OF C). 

IF U AND S2 ARE THE MEAN AN7 VARIANCE OF 

C.. . . .POLYNOMIAL IN I1 WITH EXPONENTS N,N-2,1\-4,. . . AWD IN S? WITH 

DIMENSION C(l) 

NTERMS=IFIX (FLOAT(N)/2.0+1.001) 

IF (NTERMS-1) 600,600,lOn 

inn DO 200 K=Z,NTERMS 

300 C(K)=O.Q 

DO 500 NN=2,N 

SAVEl=l.O 

DO 400 K=;Z,NTERMS 

L=NN+3-K*2 

IF (L) v-m,5noYmo 

300 IF (SAVEl-.lEt30) 350,500,500 

C(K)=SAVE2+SAVEl*FLOAT(L) 

350 SAVE;Z=C(KJ 

400 SAVEl=SAVE? 

500 CONTINUE 

mn c(i)=i.c) 

R ETUR N 

END 

SURROUT I WE CLPOP (X , V, NR ,N, IT, ICONF , ENO, ENT, W,C , CC , CLL , CLll , POPL, 
27 

fr POPIJ,PR) 

C.....RETURNS UPPER ANT) LOWER CONFIDENCE LIMITS CLll AND CLL ON THE 

c FACTOR OF INCREASE AND POPU AND POPL ON THE TOTAL POPULATION 

C EFlT FOR PROGRAM SPP. THE PROBARILITY PR THAT THE FACTOR OF 

C.. . . .REQUIRED SURROUTINES: M(IMENT,DU,DS2,PMPY,PSIIB,PILD,PQSD,CIlMPR,NRINV 

C.....COMPUTE TOTAL POPULATION EXPECTATION AND VARIANCE 

c INCREASE IS GREATER THAN 1.n IS ALSO COMPUTED. 

DIMENSJON X(l),V(NR,l),W(l),C(l) ,CC(l) 

ENT=O. 0 

VAR=fl. O 

ENT=ENT+X ( I ) 

DO lOOn J=l,N 

no ioon I=~,N 

inno VAR=VAR+V ( I , J) 

IF (VAR) inzn,iozn,iono 

1020 PR=-100.O 

RETURN 

FPR=ENT/ENO 

1 n40 VAR=VAR/ ( E NO*E NO) 

IF (FPR) 1020,1020~10m 

C.....COMPUTE COEFFICIENTS OF MOMENT POLYNOMIALS 

iom IT?=IT*~ 

CALL MOMENT (IT,C, InIM1) 

CALL PMPY (W, IDIMW,C, IDIMl ,Cy IDIM1) 

CALL MOMENT ( IT? , CC , ID IM2) 

CALL PSUR (CC,IDIM?,CC,IDIM2,W,IDIMW) 

FT=FLOAT( IT ) 

R=FPR**( 1. n/FT) 

D=FT*FPR/R 

S2=VAR/ (D*D) 

T=FPR-(FT-l.n)*FT/2.O*S2*R**(IT-2) 

C.....INITIAL GUESS AT ROOTS BY TRUNCATED TAYLOR SERIES 

IF (TI m,m,65 

60 11=3 .O 

65 Il=T**( l.O/FT) 

GO TO 70 

C.....IMPROVE ESTIMATES OF ROOTS BY ITERATION 

C.....SEARCH JS TERMINATED WHEN PERCENTAGE CHANGE BETWEEN SUCCESSIVE 

C.....ESTIMATIONS IS LESS THAN EPS FOR ROTH VARIABLES. 

70 E ps=. nnooni 

IFLAG=O 

UST=U 

S7ST=S2 

C.. . . .FIRST FUNCTION 

80 no 3nn K=~,?o 

ICUT=I 0 IM1 

CALL DU ( I T , I C UT , 11 , S7 , C , W , F PR , FU , DF DU ) 

IF (FIJI 9n,200,90 

90 IF (PFI)U) 120,100,120 

100 SLOPEl=O.fl 

B1=S2 

GO TO 200 

120 D=FU/DFDU 

Ul=ll-P 
28 

IF (PI 130,140,130 

IF (IJFDSZ) 160,140,160 

130 CALL DS2 (IT, ICIJT, IDIMl ,U,SZ,C,W,FPR ,FS2,DFDS2) 

140 SLOPEl=ln.**lO 

cI0 TO 180 

160 Dn=FS2/DFDS7 

SLOPEl=DD/( -D) 

180 Bl=S2-SLOPEl *Ul 

C.. . . .SECOND FUNCTION 

200 ICIlT=IDIM2 

CALL DU ( I T2, ICUT ,U, S7 ,CC , W ,VAR ,GU,@GDU) 

IF (GU) 220,2m,220 

270 IF (DGDII) 240,230,240 

230 SLOPE7=0.0 

B2=S2 

GO TO 290 

740 D=GU/DGDU 

112=U-n 

IF (D) 245,250,245 

745 CALL DS2 ( IT2, ICUT, IDIM7 ,\I, S2 ,CC, W ,VAR ,GS7 ,DGnSZ) 

IF (DG~SZ) 260,25n,nin 

250 sLoPEz=in.**in 

Gn TO 370 

760 D@=GSZ/DGDS2 

370 B2=S7-SLOPE?*\l7 

S L.OPE?=DD/ ( -D ) 

C.....SOLVE FOR NEW VALUES AND TEST FOR CONVERGENCE 

290 DIFF=SLOPEl-SLOPE2 

IF (DIFF) 295,310,295 

295 U=(P2-R1)/DIFF 

IF (U) 296,310,296 

296 DELTU=ARS((U-ULAST)/U) 

S2=SLOPE 2*lJ+B 2 

IF (S2) 297,310,297 

797 DELTSZ=ARS( (SZ-S2LAST)/S2) 

IF (DELTU.LT.EPS .AND. DELTSZ.LT.EPS) Gn TO 320 

UL A ST=U 

S?LAST=S2 

300 C ONT I NUE 

C.....FAILURE TI1 CONVERGE 

?in IF (IFLAG) 315,315,102n 

315 tI=lIST-. nl*UST 

S?=S2ST+.Ol*S2ST 

IFLAC=l 

GO TO 80 

C.....COMPUTE CONFIDENCE LIMITS 

320 IF (s2) io2o,inm,340 

CONF=FLOAT( ICONF) /inn. 

340 STD=SORT(S2) 

P=(l.O-CONF)/Z.n+CONF 

CALL NRINV (P,Z) 

CLU=U+Z*STD 

POPU=ENO*CLU**IT 

C L L=U - Z* S TI' 
29 

POPL=ENO*CLL**IT 

CALL CUMPR (l.O,U,STD,PR) 

R E TI JR F.l 

END 

PR=1. n-PR 

SURRCIUT I NE DIJ ( IT, JCIIT, IJ, S7 ,C, W ,CONST, FU, DFDlJ) 

C.....RETllRNS THE VALUE OF THE FUNCTION FU AND THE PARTIAL DERIVATIVE 

C.....DFDU EVALUATED AT IJ AND S7 FOR THE MOMENT POLYNOMJAL, WHOSE 

C.....COEFFICIENTS ARE GIVEN IF! VECTOR C. W IS A WORK VECTOR. 

C.....TO AVOID ARTTHMETIC OVERFLOW, HIGHER nRKR TERMS IN S2 ARE 

C.. . ..TRUNCATED WHEN THEY BECOME VERY SMALL (ICUT IS CUTOFF INDEX). 

DIMENSION C( 1) ,W( I) 

ID IM= I T+1 

J= ID IM+? 

T=l.O 

I =n 

ion I=I+I 

J=J-2 

K=J+1 

W (K)=O. 0 

700 W(J)=C( I)*T 

T=T*S7 

IF (J) 3n0,3n0,7n0 

IF (1-1) 220,2m,m 

740 IF (N(J)-F) ~~,xo,I~o 

280 w(J)=n.n 

720 F=W(J)/lOOOO. 

GO TO lclO 

260 ICIIT=I 

J=J - 1 

J=J - 1 

IF (J) mn,m,m 

300 W(l)=W(l)-CONST 

CALL PI LD (FU, DFDU , [I, W, IDIM! 

RETURN 

END 

SURROUTINE DS? (IT,~CIlT,I~IM,U,S7,C,W,C~~ST,FS?,~FDS2) 

C.....RETURNS THE VALUE OF THE FUNCTION FS2 AND THE PARTIAL DERIVATIVE 

C.....DFDS7 EVALUATED AT IJ AND S7 FOR THE MOMENT POLYNOMIAL, WHOSE 

C.....CnEFFICIENTS ARE GIVEN IN VECTOR C. W IS A WORK VECTOR. 

C.....IDIM IS ONE MPRE THAFI THE DEGREE OF THE POLYNOMIAL IN S2. 

C.....TO AVOID ARITHMETIC OVERFLOW, HIGHER ORDER TERMS IN S2 ARE 

C.. ... TRUNCATED WHEN THEY BECOME VERY SMALL (ICUT IS CUTOFF INDEX). 

DIMENSION C(l),W(l) 

U 2=lJ*U 

IJIT=IJ**IT 

W( I)=C( I)*UIT 

DO 7nn I=~,ICIIT 

200 UIT=UIT/lJ2 

600 l=I+1 

mn IF (I-IDIM) 400,500,500 

w( I )=n.n 
30 

GO TO 300 

500 W(l)=W(l)-CONST 

CALL PILD (FS2 ,DF@S2, S2 ,W, IDIM) 

RETURN 

END 

SUBROUTINE ClfMPR (X, XBAR,STD,PR) 

C.....RETURNS THE AREA PR UNDER THE NORMAL PRORABILITY CURVE WITH GIVEN 

C.. . . .MEAN AND STANDARD DEVIATION LYING RELOW THE SPECIFIED VALUE X. 

C.....THIS VERSION USES HASTING'S APPROXIMATION FOR THE INTEGRAL 

C.. . . .OF A NORMAL PROBABILITY DENSITY FUNCTION. 

DATA ~/n.3989423/~/n.z3i~4ig/~i/n.3ig3~i~30/~~/-~.35~5~~787/ 

* m/i. 7tx477937/~4/-i. ~1255978/~5/1.33n274429/ 

Y=X-XBAR 

Z=Y /STD 

PR=A*EXP (-Z*Z/2. D)*T*( Bl+T*( R2tT*( R3tT*( B4tT*R5) ) ) ) 

T=l. n/( l.O+P*ABS( Z) ) 

IF (Y) IO~O,I~OO,~OO 

9nn PR=l.O-PR 

ion0 RETURN 

END 

SURROUTINE NRINV (P,Z) 

C.. . . . INVERSE OF CUMULATIVE DISTRIBUTION FUNCTION FOR GAUSSIAN WITH 

C ZERO MEAN AND UNIT VARIANCE. USES HASTING APROXIMATION TO 

C FIND Z FOR THE FllNCTION Z=G(P), WHERE P=F(Z) IS THE INTEGRAL 

C OF THE NORMAL CURVE FROM MINUS INFINITY TO 2. 

C.....THE TRANSFORMATION AT STATEMENT 1070 AVOIDS AN ILLEGAL ARCUMENT 

C CONDITION (LOG OF ZERO) IN THE FOLLOWING STATEMENT FOR VALUES 

DATA AZ/2.515517/A1/~.~~2853/A2/n.OlO328/Rl/l.~327~~/~?/~.1~9~~9/ 

c.....oF P NEAR ZERO OR ONE (SPECIFIC FOR IMS !jnnn). 

ff ~3/n.nni308/ 

PM5=P-O. 5 

IF (PMF;) inm,in50,i040 

in40 Q=~.o-P 

io50 z=o.n 

in70 O= (otn. nnonooi ) /I. noonno1 

GO TO 1070 

GO TO in90 

1060 O=P 

T =SOR T ( - 2. fi*A LOG ( O ) ) 

T2=T*T 

T3=T2*T 

Z=T-(AZ+A1*T+A2*T2)/(1.O+RlfT+R2*T2+R3*T3) 

IF (PM~) in80,1050,in9n 

imn z=-z 

in90 RETURN 

END 

SUR ROUT I NE PMP Y ( Z , In IMZ , X , I D IMX , Y , I D I MY ) 

C.....MlJLTIPLIES POLYNOMIAL X RY POLYNOMIAL Y. RESULTANT VECTOR Z 

C.....CANMOT RE IN THE SAME LOCATION AS EITHER X OR Y. DIMENSIOF'S 

C.....OF VECTORS OF COEFFICIENTS ARE ONE MORE THAN THE DEGREE OF THE 

C.....POLYNOMIAL. COEFFICIENTS ARE ORDERED FROM SMALLEST TO LARGEST. 
31 

C.....THIS VERSION IS MODIFIED TO AVOID ARITHMETIC OVERFLOW WITH 

C.....VFRY LARGE COEFFICIENTS COMPUTED WITH SUBROUTINE MOMENT. 

DIMENSION Z( 1) ,X( 1) ,Y (1) 

IF (IDIMX*IDIMY) in,in,~n in IDIMZ=O 

GO TO 50 

2n IDIMZ=IDIMX+IPIMY-1 

KMAX= I D IMZ 

DO 30 I=l,IDIMZ 

DO 32 I=l,IDIMX 

30 z(I)=n.n 

IF (X(I)-.lE+lR) 3?,32,34 

32 CONTINUE 

GO TO 35 

34 KMAX=?*I-2 

35 DO 40 I=l,IDIMX 

DO 40 J=l,IDIMY 

IF (K-KMAX) 38,38,36 

K=I +J - 1 

36 Z(K)=.lE+20 

GO TO 40 

40 CONTINUE 

38 Z(K)=Z(K)+X( I)*Y(J) 

50 RETURN 

END 

SUBROUTINE PSUR ( Z, IDIMZ ,X , IDIMX, Y, IPIMY) 

C.....SllBTRACTS POLYNOMIAL Y FROM POLYNOMIAL X. RESULTANT VECTOR Z 

C.....MAY RE IN THE SAME LOCATION AS EITHER VECTOR X OR Y ONLY IF 

C.....THE DIMENSION OF THAT VECTOR IS NOT LESS THAN THE OTHER INPUT 

C.....VECTOR. COEFFICIENTS ARE ORDERED FROM SMALLEST TO LARGEST. 

C.....PIMENSIONS OF X AND Y ARE ONE MORE THAN THE DEGREE OF THE 

C.....POLYNnMIAL. 

DIMENSION Z(l) ,X(1) ,Y(1) 

NDIM=IDIMX 

IF (IDIMX-IDIMY) 10,2n,20 

in NDIM=IDIMY 

2n IF (NDIM) 9n,9nY3n 

IF ( I-IDIMX) 40,40,fi0 

40 IF (I-IDIMY) 50,50,70 

50 z ( I )=x ( I ) -Y ( I ) 

60 Z(I)=-Y(1) 

GO TO 8n 

70 Z(I)=X(I) 

30 DO 80 I=l,NDIM 

GO TO 80 

80 CONTINUE 

90 IPIMZ=NDIM 

RETURN 

END 

SUBROUTINE PILD (POLY,DER,ARGUM,X,IDIMX) 

C.. . . .EVALUATES POLYNOMIAL AND ITS FIRST DERIVATIVE FOR A GIVEN 

C.....ARGUMENT BY MEANS OF OUAPRATIC SYNTHETIC DIVISION (POSD). 
32 

C.. . . .VECTOR X CONTAINS THE COEFFICIENTS OF THE POLYNOMIAL. 

DIMENSION X(l) 

P=ARGUM+ARGUM 

CALL POSD (DER ,POLY ,P ,Q, X, IDIMX) 

POLY=ARGUM*PER+POLY 

RETURN 

END 

O=- AR C,IJM*AR GUM 

SUBROUTINE POSD (A,B,P,Q,X,IDIMX) 

DIMENSION X(l) 

A=O. 

B=O. 

J= I D IMX 

1 IF (J) 3,3,2 

2 Z=P*A+B 

R=O*A+X (J ) 

A=Z 

GO TO 1 

3 RETURN 

END 

C.....OUAPRATIC SYNTHETIC DIVISION 

J=J-1 
33 

APPENDIX 4. Listing of program SLT. 

PRPGRAM SLT -- STOCHASTIC LIFE TARLE SIMULATOR 

C 

C**** POPULATION PRPJEI:TIQF1 USING STOCHASTIC LESLIE MATRICES *********** 

C INPUT PATA AFJD CONTROL PARAMETERS ARE SPECIFIED ON INPUT FILE 

C '1NPUT.SLT'. FOR INSTRUCTIONS, REFER TO FILE SLTGUIDE. 

C.....THE MEAFJ (DETERMINISTIC) VALUES FPR PROJECTION ARE GIVEN RY THE 

C 

C VERSIONS ADD RANDOMNESS IN ACCORDANCE WITH VARIANCES AND COVARIANCES 

C SPECIFIFn IN THE COVARIANCE MATRIX S. 

C.....SPECIAL OPTIONS ARE: 

I: IOPl=P INITIAL POPULATION VECTOR (RNX) IS SPECIFIED IN DATA 

C IOP1>0 INITIAL POPULATION VECTOR IS CALCULATED AS THE STARLE ACE 

C DISTRIPIITION; THE TOTAL INITIAL POPIJLATION MUST RE SPECIFIED 

C IOP?=P FULL COVARIANCE STRllCTURE 

C IPP2>P SIMPLIFIED COVARIANCE STRUCTURE; A SINGLF COVARIANCE RETWEFN 

C ALL FECUNDITIES ANI) ALL SURVIVAL RATES IS SPECIFIED 

C 

SCHEDULES OF FECIINDITY (RX) AND SURVIVAL RATES (PX). THE STOCHASTIC 

DIMENSION RNX(2O) ,ENX( 20) , XBAR( 21)) ,XCVAR( 211)) ,RX (20) ,PX( 20) , S(820) 

DIMENSION RT(2O) ,PT(2P) ,RTRAR(2O) ,PTBAR(2fl) ,RTVAR(ZP) ,PTVAR(20) 

DIMENSION VX(2P) ,VXRAR(?O) ,VXVAR(?0) ,@ERDM(?O) ,DRDMM(20) ,DRDMV(20) 

DIMENSION DERDP(2n) ,DRPPM(70) ,PRDPV(20) ,ELX(20) ,PHIX(20) 

DIMENSION BPT(4O) ,RPTRAR(dfl) ,RPCV(820) ,W1(41)) ,W2(41)) ,V(1600) 

DIMENSION ITABl(50) ,ITAR2(50) 

COMMON ZERO(dO) 

DATA MM/1)/NRAD/O/ITAB1/5O*~/ITAR~/~O*O/NO~S/fl/ZERO/4(1*O.0/ 

DATA IP/fl/ I CLL/1)/ ICLU/O/ 

CALL OPEPJ (I , ' INPUT SLT' ,P) 

CALL OPEN (6,'OIlTPIJT SLT',O) 

CALL LOAD (N, FIT, MR, IOPl , IOP?,LVN, ISEED, RNX,RX , PX, S,BPCV, KK, W I@THl , 

NRNT=NR*NT 

N 2= N* 7 

ICNT=2O/( PJ*FJT)*10(1 

IF (ICNT .EO. 0) ICNT=11) 

C.. . . .COMPUTE LIFE TABLE PARAMETERS 

CALL PTOL (PX, FLX, N) 

CALL LMPH I (F: LX ,BX, PHI X ,I!) 

CALL MMNT (PHIX,RZ ,T, S?,N) 

CALL ERNWTA (PHIX,RZ,T,S?,FR,N) 

CALL FISH (PHI X,ELX,VX,ER,N) 

CALL DERV (ELX,VX,T,ER,DERPM,DERDP,N) 

C.. . ..COMPUTE THE STABLE ACE PISTRIRIITION 

C.....SET PARAMETERS AKD OPTIONS 

P RL1, W IDTH? ,BL? ,RL, RII) 

IF (10~1) ii4o,ii4o,iinn 

11nO P=Rt'JX(l) 

CALL STABL (ER,ELX,RNX,N) 

no 1120 l=l,N 

RNX( I)=RNX( I)*P 

1120 COWTINUE 

C.. . . .INITIALIZE POPIILATION VECTOR AND PROJECT 

iim rm 3200 I= 1,~' 
34 

ENX( I)=RNX( I) 

DO 1210 NTT=l,NT 

CALL PRJCT (ENX,PX,BX,N) 

CALL EREAL (RNX , ENX , N,NT, FI , ERT) 

CALL OUT (1 , N , RNX , ENX , ZERO, RX , ZERO, PX , ZEliO, VX , ZERO, DERDM ,DERDP , 

1700 CONTINUE 

1210 CONTINUE 

C.. ... OUTPUT OF DETERMINISTIC PROJECTION AND ECHO OF JNPUT 

* S,RPCV,ER,O.,RZ,0.,T,O.,S2,O.,F1,O.,ERT,O.,1P,NR,NT,1,O, * LIJN, JOPl ,KK, W IDTH1 ,BL.l ,WIDTH2 , BL?, ITARl , ITAB? , RL, RU, 0,O) 

C.. . . .EIGENANALYSIS OF COVARIANCE MATRIX (REVERSE PCA) 

C.. . . .MESSAGE TO OPERATOR 

WRITE (LUN,13) 

CALL EIGEN (s,v,w,n) 

13 FORMAT (//SX,'BEGINNING STOCHASTIC PROJECTIONS NOW . . .I/) 

C.. . . .STOCHASTIC PROJECTIONS 

NDX=O 

PO 3000 NRR=l,NR 

C.....INITIALIZE POPULATION VECTOR 

PO i8no I=I,N 

ENX( I)=RNX( I) 

1800 CONTINUE 

C.....GENERATE STOCHASTIC VECTORS, ANALYZE AN@ PROJECT 

DO 1840 NTT=l,NT 

CALL XRGEN (RPT,S,V,WI ,W2,N7,ISEED) 

CALL RPOUT (BPT,RT,PT,N) 

CALL PTOL (PT,FLX,N) 

CALL LMPH I (ELX,BT, PHI X , N) 

CALL MMNT (PHIX,RZ,T,S2,N) 

CALL ERNWTA (PHI X ,RZ ,T, S2, ER ,P) 

IF (ER+QO.O) 1825,1825,183O 

1875 NBAP=NRAn+l 

GO TO 1820 

1830 CALL FISH (PHI X , ELX,VX , ER ,FJ) 

CALL DERV (ELX,VX,T,ER,DER@M,I?ERDP,N) 

CALL PRJC T (E NX , PT , RT , h' ) 

NDX=NDX+l 

CALL MVAC (ER,ERBAR,ERVAR,NnX,NRNT) 

CALL MVAC (T,TRAR , TVAR , NDX, NRNT) 

CALL MVAC (S? , S2RAR , S7VAR , NDX , NRNT) 

CALL MVAC (RZ ,RZRAR,RZVAR,NDX,NRNT) 

CALL MVVAC (VX,VXBAR,VXVAR,N,NDX,NRNT) 

CALL MV VAC ( n E R DM , RR PMM , RR DMV , N , PIn X , FIR NT ) 

CALL MVVAC (PERDP,PRPPM,DRPPV,N,NDX,RRNT) 

CALL CVAC (RPT,BPTRAR,RPCV,N2,NDXy~lRNT) 

CALL CVAC (ENX,XBAR,XCVAR,N,h'RR,NR) 

CALL EREAL (RNX , ENX ,N , NT , F I , ERT) 

CALL MVAC (F J , FIRAR , FIVAR ,NRR,NR) 

CALL MVAC (ERT,ERTRR,ERTVR,NRR," 

CALL TALLY (F I , KK , WI@THl , RLl , NORS, ITARl ) 

CALL TALLY (ERT, KK , WIDTH2, BL?, NORS, ITAR2) 

1820 CALL BPIN (BX,PX,BPT,N) 

C.....ACCUMULATE SUMS FOR MEANS AND VARIANCES 

18411 CONTINUE 
35 

CALL COWT (NRR,ICNT,MM,LUN) 

IF (ERT-1.0) 187o,1870,186~ 

1860 IP=IP+l 

1880 I C I.L= I C L L t 1 

1900 I c LU= I c LlI+1 

1870 IF (ERT-RL) imoyi8m,iR9n 

mn IF (ERT-RII) mnc,im,im 

mn CONTINUE 

C. .OUTPIIT OF STOCHASTIC PRflJECTI0NS 

CALL RPOUT (RPTBAR , RTRAR , PTBAR, P!) 

CALL VAROUT (RPCV,RTVAR,PTVAR,B) 

CALL OIJT (?,N,RWX,XRAR,XCVAR,BTRAR,RTVAR ,PTBAR, TVAR,VXBAR,VXVAR, 

* DR~MM,~RDPM,RPCV,S,ER~AR,ERVAR,RZRARyRZVARyTBARyTVARyS2BARy 

* S2VAR,FIRAR,FIVAR,ERT~RyERTVR,IP,NR,NTyN~Xy~~BADyL~JNyIOPl,KKy 

* WI@THl,RLI,WI~TH2,BL7,ITABl ,ITAR2,RL,RII,ICLL,ICLU) 

C. . .DETERMINISTIC PROJECTION WITH ACHIEVE@ MEAN VALUES 

C.. . . .COMPUTE LIFE TABLE PARAMETERS FOR ACHIEVED MEAN VALUES 

CALL PTOL (PTBAR,ELX,W) 

CALL LMPHI (ELX,BTRAR ,PHIX,FI) 

CALL MMNT (PHIX,RZ ,T, S2,N) 

CALL ERNWTA (PHIX,RZ,T,S?,FR,N) 

CALL FISH (PHIX,ELX,VX ,ER,N) 

CALL nERV (E LX, VX , T , ER ,DEREM ,DER@P N) 

C.. . . .COMPUTE WEW STABLE AGE DISTR IRUTION 

IF (10~1) 2080,2nm,7040 

mm CALL STABL (ER,ELX,RNX,N) 

no 2nm I=I,N 

2nm CPFITINUE 

2080 no ?ion I=~,N 

7100 CONTINUE 

RNX( I)=RNX( I)*P 

C.. ... INITIALIZE POPULATION VECTOR ANI, PROJECT 

ENX( I)=RNX( I) 

DO 7120 NTT=l,NT 

CALL PRJCT (ENX,PTBAR , RTBAR ,N) 

CALL EREAL (RNX,ENX,N,NT,FI, ERT) 

CALL OUT (3,N,RNX,ENX,ZERO,RTBAR,ZERO,PTBAR,ZEROyVXyZER~yDERDMy 

* @ER~P,S,RPCV,ER,~.,RZ,O.,T,O.,S2,O.,F1,O.,ERT,O.,1P,1,NT,1, * 

21 20 CONTINUE 

C.. . . .OUTPUT OF SECOND DETERMINISTIC PROJECTION 

0, LUN , n,KK, W InTHl , RLJ , W IDTHZ,BL2, ITABl , ITAR2 , RL, RU, O,0) 

CALL EXIT 

EN@ 

SUR ROUT I NE OIJT ( I CODE , N , RNX , XRAR , XC VAR , RXR AR , RXVAR , PXB AR , P XVA R , * VXRAR,VXVAR,DR~MM,~RDPM,C~VAR,W,ER,ERVARyRZyRZVARyTyTVARyS2, 

* S2VAR,F1,F1VAR,ERT,ERTVAR,1P,NR,~JT,NDX,FIRA~,LUN,1~P1,KK, 

* W1~TH1,BL1,W1~TH2,BL2,1TAR1,1TAR2,RL,RU, ICLL, ICLU) 

C*****OUTPUT FOR STOCHASTIC LIFE TABLE SIMULATION OF NT TIME STEPS WITH 

C NR REPETITIONS. THE VECTORS ARE OF LENGTI; V WITH MEAFIS CALCULATED 

C OVER NDX VALUES. W IS A WORK VECTOR DIMEVSIONED .GE. COVAR. 

DIMENSION RNX(1) ,XRAR(l),XCVAR(l) ,BXRAR(l) ,RXVAR(l) ,PXRAR(l) ,W(1) 

nIMENSION PXVAR( 1) ,VXBAR( 1) ,VXVAR( 1) ,PRDMM( 1) ,DRDPM( 1) ,COVAR( 1) 

DIMENSION ITABl(l),ITAB2(1) 
36 

C.....HEAPING FOR VITAL RATES 

Gn TO (iioo,~onn,iiw) ,mDE 

WRITE (6,in) 

1100 WRITE (LUN,lO) 

10 FORMAT (1H1///9XY'* * * * * STOCHASTIC LIFE TABLE SIMULATION * * * 

* * *'///5Xy'l. INPUT'//5XY'THE FOLLOWING VITAL PARAMETERS HAVE RE 

*EN SPEC IF IED: '/) 

GO TO 1160 

1150 WRITE (LUN,l5) 

WRITE (6,15) 

15 FORMAT (1H1///5X , '4. DETERMINISTIC POPULATION PROJECT10N1///5X , *'THE ACHIEVEP MEAPil VITAL RATES PRODUCE THE FOLLOWING ASSOCIATED' 

*lX,'PARAMETERS:'/) 

C.. . . .TARLE OF VITAL RATES 

1160 WRITE (LUNJO) 

WRITE (6,3O) 

30 FORMAT (/6X , 'AGE ' ,4X , 'BIRTH' ,5XY 'SURVIVAL' ,7X, 'REPROPUCTIVE '/Fix, 

*'CLASS' ,4X, ' RATE' , 7X , ' RATE' , 7X , 'VALUE' ,5X, 'DER/PM' ,5X, 'DER/PP'/) 

DO 1170 I=I,N 

WRITE (LUN,?5) I ,RXRAR( I) ,PXBAR( I) ,VXRAR( I) ,[IRDMM( I) ,DRPPM( I) 

WRITE (6,35) I ,RXRAR( I) ,PXRAR( I) ,VXRAR( I) ,DRDMM( I) ,PRDPM( I) 

35 FPRMAT (6XYI2,5F11.5) 

117n CONTINUE 

WRITE (LIJN , 40) ER, RZ , T , S2 

WRITE (6,4O) ER,RZ,T,S? 

40 FPRMAT (/5X,'ASYMPTOTIC FACTOR OF I~CREASE',lX,FlO.4/ 

*5X,'NET REPLACEMENT RATE',lOXYF1O.4/ 

*5X,'GENERATION TIME',lSX,Fl0.4/ 

*5X,'VAR. IN AGE OF REPRODUCTION' ,3X,FlO.d) 

C.. . . .POPULATION VECTOR TABLE 

Gn Tn (i200,~000,i2m) ,IC~DE 

'1201) WRITE (LIlN,45) NT 

WRITE (6,45) NT 

45 FORMAT (///5X,'USING THESE VITAL RATES, THE INITIAL POPULATION IS 

* PROJECTED FOR',I3,' TIME STEPS AS FOLLOWS:') 

IF (10~1) i3noYi3o~,i2io 

no 1220 I=~,N 

inn s=n.n 

S=StRNX( I) 

WRITE (LUN,46) S 

WRITE (6,46) S 

1320 CONTINUE 

46 FORMAT (9X,' (YOU HAVE SPECIFIED THAT THE INITIAL POPULATION SHOULD 

* RE IN STABLE AGE ~ISTRIBIITION'/9X,'WITH A TOTAL OF',F10.11, 

*' INDIVIDUALS) I) 

Go TO 1300 

1260 WRITE (LUN,~~) NT 

WRITE (6,48) NT 

48 FORMAT (///5X,'IJSING THESE ACHIEVED MEAN VITAL RATES, THE INITIAL 

* PnPIILATION'/SX,'IS PROJECTED FOR',I3,' TIME STEPS AS FOLLOWS:') 

i3nn WRITE (LIJNJO) 

WRITE (6,5O) 

50 FORMAT (//l~X,'A~E',~X,'I~lITIAL', 9X,'FINAL'/ 

*9X,'CLASS',3X,'POPULAT10N',5X,'POP1lLAT1ON'/) 
37 

no 1310 I=l,N 

WR ITE (LUN, 57) I, RNX( I ) , XRAR( I ) 

WRITE (6,FiZ) I,RNX(I),XRAR(I) 

52 FORMAT (lOX,I?,?F15.2) 

1310 CONTINUE 

WRITE (LUFI,53) ERT,FI 

WRITE (6,53) ERT,FI 

53 FORMAT (/SX,'FINAL/INITIAL POPULATION',6X,FlO.4/ 

*SX,'REALIZED FACTOR OF INCREASE' ,3X,FlO.4) 

GO TO (I~OP,?OOO,~OOO) ,ICODE 

C.....LIST SPECIFIED COVARIANCE MATRIX 

1400 N?=N*? 

blR I TF ( LUb!, 60) Pi, N 

60 FORMAT (lHI///!iX,'?. 

WRITE (6,601 N,V 

STOCHASTlC POPULATION PROJECTION'// 

*Fix, 'THE FOLLOWING COVARIANCE MATRIX HAS BEEN SPECIFIED: '/ 

*lnX,' (FIRST' ,I3,' COLUMNS ARE RIRTH RATES, SECOND',I3, 

*' ARE SURVIVAL RATES)'/) 

CALL WTR I (N? ,COVAR ,W ,LUW, 2) 

CALL WTRI (N?,CQVAR,W,6,?) 

GO TO 3000 

C.....STOCHASTIC PROJECTION OUTPUT 

C.. ... POPULATION VECTOR TABLE 

WRITE (6,7O) NT,NR 

znnn WRITE (LUN,~~) NT,NR 

70 FORMAT (///5X,'llSING THE MEAN VITAL RATES WITH THIS COVARIANCE' 

*lX,'STRIJCTURE,'/SX,'THE INITIAL POPULATION IS PROJECTED FOR',I3, 

*' TIME STEPS',I6,' TIMES.') 

WRITE (LUN,72) 

WRITE (6,72) 

72 FORMAT (/lPX,'AGE',6X,'INITIAL',FiX,'(MEAM) FINAL'/ 

*9X, 'CLASS' ,3X, ' POPULATION' ,5X, 'POPULATION'/) 

no 7220 J=l,N 

WRITE (LUN, 52) I ,RNX( I ) ,XBAR( I) 

WRITE (6,5?) I ,RNX( I) ,XBAR( I) 

mn CONTINUE 

C.. ... POPULATION VECTOR COVARIANCE MATRIX 

WRITE (LUN,74) 

WRITE (6,74) 

CALL WTRI (N,XCVAR,W,LUN,l) 

CALL WTRI (N,XCVAR,W,6,1) 

WRITE (LllN,75) ERT,ERTVAR,FI,FIVAR 

WRITE (6,75) ERT, ERTVAR, F I, FIVAR 

74 FORMPT (//SX,'COVARIAMCE MATRIX FOR FINAL POPULATION VECTOR'/) 

75 FORMAT (/40x, 'MEAN' ,~X,'VARIANCE'// 

*5X,'FINAL/INITIAL POPULATION',6X,FlO.4,FlO.5/ 

*5X, 'REALIZED FACTOR OF INCREASE' ,3X,Fl0.4,FlO. 5) 

WRITE (LlJN,77) 

WRITE (6,77) 

77 FORMAT (//5X,'FRENENCY DISTRIRUTInN FOR FACTOR OF INCREASE'// 

*57X,'REALIZED'/14X,'LOWER ROIIND',?X,'FINAL/INITIAL', 3X,'LOWER' 

*1X, 'ROUND' ,2X, 'FACTOR OF ' /ax, ' INTERVAL' , 2X, 'OF INTERVAL' ,3X, 

*'POPULATION', 5X, 'OF INTERVAL' ,3X, ' INCREASE'/) 

Rl=RLl 
38 

R2=RL7 

DO 7300 I=l,KK 

WR ITE (LUN, 78) I , B2 , ITAB2( I ) ,R1 , ITABl ( I ) 

WRITE (6 , 78) I ,B2 , ITA62( I ) , 61, ITABl ( I ) 

FORMAT (5X , I7 , 5X , ' <' ,F8.4,110, 1OX , ' <' ,F8.4,IlO) 

GO TO 7300 

WRITE (LUN, 79) I ,B2, ITAR2( I ) ,B1 , ITARl( I ) 

WRITE (6,79) I ,62 , ITA62( I) ,61 , ITA81 (I) 

FORMAT (5X , I7 ,Fix ,F8.4 , IlC),llX , F8.4,IlO) 

Rl=Bl+WIDTHl 

R2=82+W IDTH2 

C ONT I NUE 

FNR=F LOAT( NR ) 

PRl=FLOAT( IP)/FNR 

PR 2=F LOAT( IC LL) /FNR 

PR 3=F LOAT( IC LU ) /FNR 

WRITE (LUN,76) PR1,RLYPR2,RU,PR3 

IF (1-1) 2240,224nY2260 

WRITE (6 , 76) PR1 , RL , PR2 , RU, PR3 

FORMAT (/SX,'PROPORTION OF PROJECTInNS WITH FACTOR nF INCREASE', 

*' GREATER THAN 1.0 IS',F7.4/5XY'PROPORTION OF FINAL/INITIAL', 

*' POPIJLATION RATIOS RELOW',FR.4,' IS',F7.4/5XY'PROPORTION OF', 

*I FINAL/INITIAL POPULATION RATInS ABOVE',F8.2,' IS',F7.4) 

C.....HEADING FOR PART 3 - STATISTICAL SUMMARY 

WRITE (LUN,80) NDX,NT,NR 

WRITE (6 ,8n) NDX,NT,NR 

an FORMAT (IHI///~X,'~. STATISTICAL PROFILES OF DEMOGRAPHIC PARAMETE 

*RS USED IN STOCHASTIC POPULATION PROJECTIONS'// 

*lflX,'THESE STATISTICS ARE COMPUTER ON THE BASIS OF' ,16,' VALUES'/ 

*lOX,'(',J?,' TIME STEPS,' ,15,' REPETITIONS)') 

WRITE (LIlN,81) 

WRITE (6,81) 

81 FORMAT (//6X,'A~E',lOX,'BIRTH RATE',RX,'SURVIVAL RATE',4XY'REPRODll 

*CTIVE VALUE',4X,'~ER/PM',~Xy'DER/DP'/ 

*5X , 'CLASS' ,3( 7X , 'MEAN' ,6X , 'VAR' ) ,7X , 'MEAN' ,5X, 'MEAW' /) 

DO 2310 I=,W 

WRITE (LUN,82) I,BXBAR(I),BXVAR(I),PXBAR(I),PXVAR(I),VXRAR(I) , 

WR I TE (6 , 82) I , BXR AR ( I ) , BXVA R ( I ) , P XRAR ( I ) , P XVAR ( I ) , VXR AR ( I ) , 

* 

* 

VXVAR(1) ,PRDMM( I) ,DRPPM( I) 

VXVAR( I) ,nRPMM( I) ,DRDPM( I) 

82 F ORMAT (6 X , I 2 , 3X , 4 ( F 1 I. 5 , F9.5) ) 

2310 CONTINUE 

C.....WRITE OTHER STATISTICS 

WRITE (LUN,85) ER,ERVAR,RZ,RZVAR,T,TVAR,S2,S3VAR 

WRITE (6,85) ER,ERVAR,RZ,RZVAR,T,TVAR,S7,S2VAR 

85 FORMAT (/4nX,'MEAN',4X,'VAR1ANCEt// 

*5X,'ASYMPTOTIC FACTOR OF INCREASE',lX,Flfl.4,FlO.5/ 

*5X,'NET REPLACEMENT RATE',lnX,FlO.n,FlO.5/ 

*5X,'GENERATION TIME',15X,F10.4,F10.5/ 

*5X,'VAR. IN ACE OF REPRODUCTIflN',3X,Flfl.4yF~~.5) 

IF (NRAD) z400,~~00,ni~ 

2315 WRITE (LUN,Rfi) NRAll 

WRITE (6,86) NRAP 

86 FORMAT (/lOX,'***** WARNING: FACTPR OF INCREASE COULD NOT RE COMP 
*IJTFDl ,15,' TIMES') 

C.....CALCULATED COVARIANCE MATRIX FOR VITAL RATES 

2400 WRITE (LUN,87) 

WRITE (6,87) 

CALL WTR I (N?,COVAR,W ,LUN,?) 

CALL WTR I (N2 ,C0VAR ,W ,6,2) 

END 

87 FORMAT (//SX,'THE FOLLOWING COVARIANCE STRUCTURE WAS ACHIEVED: I/) 

3000 RETURN 

SIJRROUTINE WTRI (I'J,COVAR,W,LUN,IOP) 

C. .WRITES THE CONTENTS OF A SINGLE-SUBSCRIPTED, UPPER TRIAI'JGULAR 

C.....MATRIX COVAR IN N BY N TABLE FORM ON DEVICE LUN. W IS A WORK 

C.....VECTOR DIMENSIONED AT LEAST AS GREAT AS COVAR. IOP=1 FOR E 

c.. . . .FORMAT, IOP=~ FOR F FORMAT (~10.5). 

DIMENSION COVPR( 1) ,W( 1) 

COMMON ZERO(40) 

WRITE (LUN,75) (J,J=l,N) 

75 FORMAT (6xYi0Iin) 

C.....REPACK MATRIX 

L=O 

no 7250 I=l,N 

L=L+l 

I J= (J *J -J ) / 2t I 

W(L)=COVAR( IJ) 

DO 7240 ,]=I ,I'J 

7740 CONTINUE 

2250 CONTINUE 

C.....LIST 

JSTP=I) 

DO 7300 I=l,N 

ZERO( I )=0.O 

JSTRT=JSTPtl 

JSTP=JSTRT+N-I 

K=I-1 

IF (10p-1) 2280,2280,~255 

2255 IF (K) 226OY2260,2770 

2260 WRITE (LUN,77) (W(J) ,J=JSTRT,JSTP) 

2270 WRITE (LUN,77) (ZERO(J) ,J=l,K),(W(J),J=JSTRT,JSTP) 

GO TO 7300 

77 FORMAT (9X,lOFlO.5) 

2280 IF (K) 2285,2285,2290 

2285 WRITE (LUN,78) (W(3) ,J=JSTRT,JSTP) 

GO TO 7300 

7290 WRITE (LUI'J,78) (ZERO(J) ,J=l,K) ,(W(J) ,J=JSTRT,JSTP) 

2300 CONTINUE 

RETURN 

END 

GO TO 2300 

78 FORMAT (9X,10E10.4) 

SIIBROUTINE EIGEFJ( A, R ,N,MV) 

C*****EICENANALYSIS OF REAL SYMMETRIC MATRIX, A, OF ORDER N, 

c RECEIVED IF4 UPPER TRIANGIJLAR SINGLE SUBSCRIPT NOTATION. 
C.....RETURNS EIGENVALUES, LARGEST FIRST, IN FIRST N ELEMENTS OF A. 

C OPTIONALLY, RETURNS EIGENVECTORS AS COLUMNS OF MATRIX R, STORED 

C IN SINGLE SUBSCRIPT NOTATION, IN ORPER CORRESPONDING TO THE 

C ORDER OF THE EIGENVALUES. 

C.....IF THE CONTROL PARAMETER MV E0UALS ONE, EIGENVECTORS ARE 

C NOT COMPUTED, AND R NEED NOT RE DIMENSIONED LARGER THAN 1. 

C.....THE INPUT MATRIX IS DESTROYED IN COMPUTATION. 

C THIS SllRROUTINE IS SLOWER THAN THE HOUSEHOLDER ALGORITHM FOR 

C LARGE MATRICES, BUT IT IS REPORTED TO RE MORE ACClJRATE WHEN THERE 

C 

C.. . . .THE UPPER TRIANGULAR ARRAY MUST BE DIMENSIONED AT LEAST N*(N+1)/2 

C IN THE MAIN PROGRAM 

C.....MODIFIED FROM SSP TO RETURN EIGENVALUES IN FIRST N ELEMENTS OF A. 

DIMENSION A(l),R(l) 

c GENERATE IPENTITY MATRIX 

SMALL E I G E N VA LU E S . 

IF(MV-1) 1O,?5,10 

10 IO=-N 

DO 20 J=l,N 

I O= I O+N 

DO 20 I=l,N 

IJ=IO+I 

R( IJ)=O.O 

IF( I-<J) 20,15,2O 

15 R(IJ)=I.O 

20 C OFITI NIJE 

25 ANORM=I).O 

c CPMPUTE INITIAL AN@ FINAL NORMS (ANORM AND AWORMX) 

DO 35 I=l,N 

DO 35 J=I ,N 

IF( I-J) 30,35,31) 

30 I A= I+( J*J -J ) /? 

ANORM=ANORM+A( I A) *A( I A) 

35 CONTIVUE 

IF ( ANORM) 165,165,40 

40 ANORM=1.414*SORT(ANORM) 

ANRMX=ANORM*O . OOOo0l /FLOAT( N ) 

IND=O 

THR=ANORM 

C INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR 

45 THR=THR/FLOAT( N) 

50 L=l 

55 M=L+1 

60 MQ= (M*M-M) /2 

C COMPUTE SIN AND COS 

L @= (L*L -L ) /2 

LM=LtMQ 

62 IF ( ABS( A( LM) )-THR) 130,65,65 

65 IND=1 

LL=LtLO 

MM=M+MQ 

X=O. 5* (A( LL) -A( MM) ) 

68 Y=-A( LM) / SORT( A( LM)*A( LM)+X*X) 

70 Y=-Y 

IF(X) 70,75,75 
41 

75 S INX=Y/ SORT( 2.0*( 1.0+( SORT( 1 .O-Y*Y) ) ) ) 

SINX2=SINX*SINX 

COSX2=COSX*CPSX 

SINCS =SINX*COSX 

78 CQSX= SORT(1.0-SINX2) 

C ROTATE L AND M COLUMNS 

I L +N* ( L - 1 ) 

IMO=N* (M-1 ) 

DO 125 I=l,N 

IF( I-L) 8O,115,80 

80 IF( I-M) 85,115,90 

85 IM=I+MO 

90 IM=M+IQ 

IO= ( I *I -I )/? 

GO TO 95 

95 IF(I-L) io~,i0~,io5 

inn IL=I+LO 

G0 TO 110 

105 IL=L+IQ 

110 X=A( IL)*COSX-A( IM)*SINX 

A( IM)=A( I L)*S INX+A( IM) *COSX 

A( IL)=X 

115 IF(MV-1) 120,125,120 

120 ILR=ILO+I 

IMR= I MQ+ I 

X=R( ILR)*CQSX-R( IMR)*SINX 

R ( IMR)=R ( I LR)*S I NX+R ( IMR)*COSX 

R( ILR)=X 

125 CONTINUE 

X=Z.n*A(LM)*SINCS 

Y=A(LL)*COSX?tA(MM)*SINXZ-X 

X=A(LL)*SINX2+A(MM)*COSXZ+X 

A( LM)= (A( LLI-A( MM) )*S INCS+A( LM)*( COSX2-S INX2) 

A( LL)=Y 

A(MM)=X 

c TESTS FOR COMPLETION 

C TEST FOR M = LAST COLUMN 

130 IF(M-N) 135,140,135 

135 M=M+1 

GO TO F;n 

C TEST FOR L = SECOND FROM LAST COLUMN 

140 IF(L-(N-1)) 145,150,145 

145 L=L+1 

150 IF( IND-1) 160,155,160 

155 IND=n 

GO TO 55 

GO TO 50 

C COMPARE THRESHOLD WITH FINAL NORM 

c SORT EIGENVALUES AND EIGENVECTORS 

160 IF(THR-ANRMX) 165,165,45 

165 IO=-N 

00 185 I=l,N 

IO= I O+N 

LL= I +( 1*1 -I ) /2 
42 

JQ=N* ( 1-2) 

DO 185 J=I,N 

JQ=JQ+N 

IF(A(LL)-A( MM) ) 17n, 185,185 

170 X=A(LL) 

A(LL)=A(MM) 

A ( MM) =X 

MM=J+ (J*J -J ) /2 

IF(MV-1) 175,185,175 

175 DO 180 K=l,N 

I LR= I O+K 

X=R( ILR) 

R( ILR)=R( IMR) 

IMR=J Q+K 

i8n R(IMR)=X 

185 CONTINUE 

J=O 

DO 190 I=l,F! 

J=JtI 

A ( I ) =A ( J ) 

190 CONTINUE 

RETURN 

END 

C.. . . .REPACK EIGENVALUE ARRAY 

SUBROUTINE MVAC (X,XBAR,XVAR,NTT,NT) 

C*****CALCULATES THE MEAN AN@ VARIANCE OF A SCALAR X 

C BY ACCUMULATING THE SUMS AND SUMS OF SOUARES EACH TIME THE SIIRROIITINE 

C IS CALLED. THE CALL TO THIS SURROUTINE SHOIILD BE INSIDE A LOOP INDEXED 

C.....BY NTT. AFTER THE LAST PASS (NTT=NT), MEAN AND VARIANCE ARE COMPUTED. 

IF (NTT-1) 49O0,4900,5000 

4900 XBAR=O. 0 

XVAR=II. O 

5000 XRAR=XBARtX 

XVAR=XVAR+X*X 

IF (NTT-NT) !i20n,5100,5100 

sion F NT=FLOAT(NT) 

XR AR = X B A R/ F NT 

XVAR=XVAR/FNT-XBAR*XBAR 

5200 RETURN 

END 

SIIRROUTINE NRMUL (X,XBAR,STD, ISEED) 

C*****RETURNS NORMALLY DISTR IRUTED RANDOM NUMBER, X, WITH MEAN XBAR 

c AND STANPARD DEVIATION STD AT EACH CALL. 

C MULLER (1959) J. ASSOC. COMPIITING MACHIVERY 6: 376-383. 

C REF: Rnx AND MULLER (1958) AWN. MATH. STAT. 39: 610-tiii. 

DATA TP1/6.2831853072/ 

CALL RANDM (U1,ISEED) 

CALL RANDM (117, ISEED) 

X=XBAR+STD*COS (TPI*IJ?)*SORT( -2.O*ALOG( II 1) ) 

RETURN 

END 

~~i=(ui+o.nnononi)/i.~~~n~ni 
43 

SURROUTINE RANDM (X,ISEED) 

C.....RAND@M NUMBER GENERATOR. RETURNS A RANDOM NUMBER X WITH A UNIFORM 

C DISTRIRUTION RETWEEN ZERO AND UNITY, INCLUSIVE. ISEED MUST BE 

C INITIALIZED WITH A POSITIVE INTEGER IN THE RANGE 1-33767, BllT THE 

C VALUE OF ISEED MUST NOT RE RESET RETWEEN CALLS. 

C.....THI! :ENERATOR IS SPECIFIC FOR THE IMS 5000 AND IRM 1800 (6/1/81). 

COMMON MARKR,MULT,DIV 

IF (ISEED) i040,1040,iom 

1O3O MIlLT=12781 

DIV=-1.0/32767.0 

I SEED= I SEED* ( - 1 ) 

MARKR=ISEED 

1040 ISEED=ISEED*MULT-2797 

IF (ISEED) in6o,iom,mo 

1050 ISEED=ISEED-32767-1 

1060 X=FLOAT(ISEED)*DIV 

1080 MIJLT=MlILT+8 

11110 RETURN 

IF (ISEED-MARKR) llOn,1080,1100 

END 

SURROUTINE MTPRD (A,B,C,NA,MA,MB) 

B (MA RY MB) TO FORM THE MATRIX C (NA PY MR) . DIMENSION A( 1) ,R( 1) ,C( 1) 

IR=0 

IK=-MA 

DO 1050 K=l,MR 

I K= I KtMA 

DO 1040 J=l,NA 

IR= I R+l 

IB=I K 

C( IR)=fl.O 

DO 1030 I=l,MA 

JI=J I tNA 

IB=IB+l 

C( IR)=C( IR)tA(JI)*B( IB) 

C*****MATRIX MULTIPLICATION. APPLIES MATRIX A (NA RY MA) TO MATRIX 

C 

C.....THE MATRICES ARE STORED IN SINGLE SUBSCRIPT NOTATION. 

JI=J-NA 

1030 CONTINUE 

1040 CONTINUE 

1050 CONTINUE 

RETURN 

END 

SURROUTINE FISH( PHI X, ELX,VX, ER,N) 

C*****CALCIJLATES VECTOR OF REPRODUCTIVE VALUES GIVEN N-ORDER NET MATERNITY 

C.....FUNCTION, SURVIVORSHIP SCHEDULE, AND ER, THE FACTOR OF INCREASE. 

DIMENSION PHIX(1) ,ELX(l),VX(l) 

SllM=O. 0 

ERX=ER**N 

IX=Ntl 

DO in50 I=I,N 
44 

IX=IX-1 

SIJM=SUM+PHIX( IX)/ERX 

ERX=ERX/ER 

IF (ELX(1X)) 1030,1030,1040 

GO TO 1O5c) 

1030 VX(IX)=O.T) 

1040 VX( I X)=SUM*ERX/ELX( I X) 

1050 CONTINUE 

RETURN 

END 

SUBROUTINE PRJCT( ENX,PX, EMX,N) 

C*****ADVANCES A POPULATION VECTOR ONE TIME UNIT FORWARD BY APPLICATION 

C.. . . .SURVIVAL RATES. 

C OF A LESLIE MATRIX DEFINED RY N-ORDER VECTORS OF FECUNDITY AND 

DIMENSION ENX(1) ,PX(1) ,EMX(l) 

SIJM=c). O 

I X=N+l 

I XMl=N 

DO 1050 I=?,N 

I x= I x-1 

I XMl=I XM1-1 

SUM=SUM+ENX( IX)*EMX( I X) 

EMX( IX)=ENX( IXMl)*PX( IXM1) 

EFIX( l)=SUM+EMX( l)*EMX( 1) 

RETURN 

END 

1050 CONTINUE 

SIJRROUTINE PTOL( PX, ELX , N) 

DIMENSION PX(1) ,ELX(1) 

ELX( 1 )=1.O 

ELAST=1.0 

IM1=0 

I)O 1050 I=?,N 

1 PI 1 = I M1+1 

ELAST=ELAST*PX(IMl) 

ELX(I)=ELAST 

RETURN 

END 

C*****CONVERTS AN N-ORDER VECTOR OF SURVIVAL RATES TO SURVIVORSHIPS. 

in50 CONTINUE 

SIIRROUTI NE DERV( ELX , VX ,6, ER ,PERDM ,DERDP, Id) 

C*****CALCULATES THE N-ORDER VECTORS OF DERVIATIVES OF ER WITH RESPECT TO M(X) 

C ANI) P(X) FOR ALL X. 

C REOUIRES INPUT 0F THE SURVIVORSHIP VECTOR, REPRODUCTIVE VALUE VECTOR, 

C.. ... THE POPULATION GENERATION TIME, AND THE FACTOR OF INCREASE. 

DIMENSION ELX(l),VX(l),IIER~M(l),DER@P(1) 

ER lMX=ER 

IXP1=1 

DO 1050 I=l,N 

IXPl=IXPl+l 

ERlMX=ERlMX/ER 
45 

ELERG=ERlMX*ELX( I )/G 

DERDM(I)=ELERG 

IF (I-N) i04n,i06n,in40 

10411 PERDP( I)=ELERG*VX( IXP1) 

1050 CONTINUE 

ln6O DERDP(N)=n.Cl 

RETURN 

END 

SUBROUTINE MVVAC (X,XRAR,XVAR,N,NTT,NT) 

C*****COMPUTES THE VECTORS OF MEANS AND VARIAFKES OF THE N-ORDER VECTOR X 

C RY ACCUMULATING THE SUMS AND SUMS OF SOUARES EACH TIME THE SUBROUTINE 

C IS CALLED. THE CALL TO THIS SUBROUTINE SHOULD BE INSIDE A LOOP INDEXED 

C.....RY NTT. AFTER THE LAST PASS (NTT=NT), MEANS ABD VARIANCES ARE COMPUTED. 

DIMENSION X(l) ,XRAR(l) ,XVAR(l) 

IF (NTT-1) 4900,490fly5O6fl 

rl.9on DO Finon I=~,N 

Finno CONTINUE 

Fimn no t;inn I=I,N 

XRAR ( I )=O. 0 

XVAR ( I ) =O . 0 

A=X(I) 

XRAR( I )=XBAR( I )+A 

XVAR(I)=XVAR(I)tA*A 

liion CONTIN~JE 

IF (NTT-NT) 53n0,52 

Dn 5300 I=~,B 

5700 FNT=FLOAT(NT) 

S=XBAR( I ) 

XRAR ( I )=S/FNT 

200 

XVAR ( I )= ( XVAR( I )-S*S/FNT) /FNT 

53nn CONTINE 

RETURN 

END 

SURROUTIPJE CVAC (X , XRAR ,C ,N, NTT,NT) 

C**** CALCllLATES THE COVARIANCE MATRIX C OF THE N-ORDER VECTOR X RY 

C ACCIJMULATING SUMS OF CROSS PRODUCTS EACH TIME SUBROUTINE IS CALLED. 

c THE CALL TO THIS SUBROUTINE SHOULD RE INSJPE A LOOP INDEXED BY NTT. 

C WHEN NTT=NT, MEANS, VARIANCES AND COVARIANCES ARE COMPUTED. 

C.. . . THE TRIANGULAR COVARIANCE MATRIX C IS STORED IN SINGLE SURSCRIPT 

C 

C.... WHERE IJ=(J*J-J)/2tI AND 1.LE.J. THE DIMENSION OF C IS N*(N+1)/2. 

NOTATION. ELEMENT (1,J) OF THE MATRIX IS ELEMENT IJ OF THE VECTOR, 

DIMENSION X(l) ,XBAR(l) ,C(1) 

IF (NTT-1) 60OO,6000,6200 

DO 6100 J=l,N 

XRAR( J)=O. 0 

NPX=NDX+l 

6noo NDX=O 

DO f;nm I=~,J 

C(NDx)=n.n 

6ORO CONTINUE 

6101) CONTINUE 

62110 NDX=O 
DO 6240 J=l,N 

A=X(J) 

XBAR( J)=XBAR( J) +A 

PO 6220 I=l,J 

NDX=NDXtl 

C ( N D X ) = C ( N D X ) t A* X ( I ) 

6220 CONTINUE 

6240 CONTINUE 

IF (NTT-NT) 7000,~3no,6300 

6300 F NT=F LOAT( NT) 

NDX=O 

DO 6340 J=l,N 

XB AR (J ) = XB AR (J ) /FNT 

DO 6320 I=l,J 

NDX=NDX+l 

C (NPX)=C( NDX)/FNT-XRAR( I )*XRAR (J ) 

6370 COFJTINUE 

6340 CONTINUE 

70nn RETURN 

END 

SIJRROUTINE STARL (ER, ELX ,CX , N) 

C*****CALCIILATES THE STABLE AGE DISTRIRUTION CX, GIVEN THE ASYMPTOTIC 

C*****FACTOR OF INCREASE AND THE SURVIVORSHIP VECTOR OF ORDER N. 

DIMENSION ELX(1) ,CX(l) 

EMRX=l .O 

EMR=l.O/ER 

EMRX=EMRX*EMR 

SUM=SIIM+ELX( I )*EMRX 

1000 CONTINUE 

EMRX=l. 0 

EMRX=EMRX*EMR 

CX( I)=ELX( I)*EMRX/SUM 

RETURN 

END 

suM=n. n 

DO innn I=~,N 

no lion I=I,N 

iino CONTINUE 

SClRROlJTINE ERNWTA (PHI X , RZ ,T, S2 , ER , N) 

C.....THIS VERSION IS MODIFIED TO ACCEPT UNIJSUAL VITAL RATES GENERATED 

C RANDOMLY IN PROGRAM SLT. 

DIMENSION PHIX(1) 

IF (s2) ion0,inoo,in2n 

Gn TO 1030 

DO 7000 J=I,ZO 

10l)O ER=RZ**(l.O/T) 

1020 ER=EXP( (T/SZ)*(RZ**( S?/T**Z) -1.0) ) 

1030 CT=o.ooni 

ERX=1 .0 

DFDR=O. O 

FR=- 1.0 

tm low I=~,FI 
47 

ERX=ERX*ER 

IF (ERx-n.mn01) znin,~nin,in40 

104n TRM=PHIX(I)/ERX 

infio CONTINUE 

FR=FRtTRM P 

DFDR=DFDR-FLOAT(I)*TRM 

DE LT=FR/DFDR 

ABFR=ABS( FR) 

IF(ARFR-CT) 1060, 1060, 1090 

IF (ABPLT-CT) 2n50,2050,1n90 

1060 ARPLT=-ABFR/DFPR 

inw ER=ER-DELT 

2onn CONTINUE 

2051) ER=ER-DELT 

mo RETIIRR 

2010 ER=-100.0 

END 

SIJRROUTINE LMPHI (ELX, EMX ,PH I X ,N) 

C*****CONVERTS AN N-ORDER VECTOR OF SURVIVORSHIPS AND A VECTOR OF 

c FECUNDITIES TO A NET MATERNITY FUNCTION. 

DIMENSION ELX(1) ,EMX(l) ,PHIX(l) 

PHIX( I)=ELX( I)*EMX( I) 

RE TCIR N 

END 

DO 1050 I=~,N 

in50 CONTINUE 

SURROUTINE MMFlT( PHI X, RZ ,T, S?, N) 

C*****CALCULATES THE NET REPLACEMENT RATE, THE MEAN T, AND THE VARIANCE 

C S2 IN AGES OF REPRODUCTION FOR A COHORT. 

DIMENSION PHIX(1) 

RZ=n. O 

XlPHI=O.O 

X?PHI=O.O 

PHI=PHI X( I ) 

FI=FLOAT( I) 

RZ=RZtPHI 

XlPHI=XlPHI+PHI*FI 

X2PHI=X?PHItPHI*FI*FI 

DO in50 I=~,N 

1050 CONTINUE 

IF (RZ) 1n60,inm,in7n 

s?=n.n 

1060 T=1.0 

GO TI) 1100 

1070 T=XlPHI/RZ 

S?=X?PHI/RZ-T*T 

iinn RETURN 

END 

SIJRR0UTINE COUNT (N,M, I ,LUN) 

C.....DISPLAYS ON THE CRT THE PROGRESS OF A LONG PROGRAM. THE CALL TO 

c SUBROUTINE COUNT IS INSIDE A LOOP INDEXED BY N, AND THE VALUE OF 

_I_ _- 
48 

C N IS DISPLAYED ON THE SCREEN EVERY M LOOPS. I IS A COUNTER WHOSE 

C VALUE MUST NOT RE RESET RETWEEN CALLS. 

I=I+1 

IF (I-M) 1080,1060y1060 * 

in6n WRITE (LuN,in) N 

imn RETURN 

in FORMAT (~X,IR,I PROJECTIONS COMPLETED.') 

I=n 

EN@ 

SURROUTINE VAROUT (VCV,BTVAR,PTVAR,N) 

DIMENSION VCV(1) ,BTVAR(l) ,PTVAR(l) 

J=N+I 

JJ=J*(J+1)/2 

RTVAR( I)=VCV( I I) 

PTVAR( I)=VCV(JJ) 

RETURN 

EKD 

DO inon I=IJ 

II=I*( I+1)/2 

inno CONTINUE 

SIJRROUTI NE TALLY (X ,K, WIDTH , BL , WOBS, I TAR) 

C*****TALLIES THE VARIABLE X INTO K-2 CLASSES OF WIDTH WIDTH, BEGINNING 

C WITH LOWER BOUNP RL. ON OUTPUT, NORS IS THE TOTAL NUMBER OF ORSERC RE SET TO ZERO BEFORE THE FIRST CALL. 

c VATIONS TALLIED INTO THE K-ORDER VECTOR ITAR. PORS AND ITAB MUST 

DIMENSION ITAB( 1) 

R=BL-WIDTH 

DO ?fWI I=l,KMl 

R=R+W IDTH 

KMl=K-1 

IF (x-R) ~~OO,~~OO,~OOO 

2nOn CONTINUE 

I=K 

3100 ITAR( I)=ITAB( I)+1 

NOR S=NOB S+1 

RETURN 

EN@ 

SUR ROUT I NE XR GE N ( X , S , V , k'l , W? , N , I SEED) 

C*****GIVEN AN N-ORDER VECTOR OF MEANS X, AN N BY N MATRIX V, THE INVERTEP 

C TRANSPOSE OF THE MATRIX OF EIGENVECTORS OF THE COVARIANCE MATRIX, 

C ANP A VECTOR S OF EIGENVALUES (VARIANCES) , THIS SUBROUTINE GENERATES 

C A RANDOM VECTOR, RETURFED AS X, WITH THE SPECIFIED MEAP, VARIANCE 

c AND COVARIANCE PROPERTIES. 

C.....CALLS SIJRROUTINES MTPRD AMI NRMUL. ISEE@ IS THE SEEP FOR THE RANDOM 

C NIIMRER GENERATnR. WPRKSPACES W1 ANP W? MUST RE DIMENSION AT LEAST N 

C IN THE MAIN PROGRAM. 

DIMENSION X(l),S(l) ,V(1) ,W1(1) ,W3(1) 

C.. . . .GENERATE RANDOM VECTOR WITH NO COVARIANCE 

innn no lion I=~,F! 

STD=SORT( ABS( S( I ) ) ) 

CALL NRMllL (Z,O.o,STD, ISEEP) 
49 

WI(I)=Z 

1 inn CONTINUE 

C.....ADD COVARIAWCE AND MEANS 

L= 1 

CALL MTPRP (V , W1, W7 , N, PJ, L ) 

DO 12nn I=I,N 

x( I)=X( I)+W2( I) 

1200 CONTINUE 

C.. . . .TRUNCATE VALUES FOR PEMnCRAPHIC APPLICATIONS 

NH=N/? 

DO 2300 K=l ,N 

XK=X (K) 

IF (XK) 2iinymn,2i2n 

2iin x(K)=n.n 

2120 IF (K-NH) z~noY7m~,7i3n 

7130 IF (XK-1.n) 2200,23nfl,214O 

7140 X(K)=l.n 

2200 CONTINUE 

RE TlJR N 

EFJD 

SIJRROUTINE LOAD (N , HT,NR , IOPl , IOP2,LUN, I SEED, X ,RX ,PX ,COVAR ,W ,KK, 

ff WIDTH1 ,RL1 ,WIDTH?,RL7,RLYRU) 

C*****READS THE CONTROL PARAMETERS APJD INPIJT VECTORS OF INITIAL POPULATION, 

C FECUNDITY, SURVIVAL AN[! COVARIANCE FOR PROGRAM SLT. 

C.....IF IOPl>O, THE TOTAL INJTIAL POPULATION SHOULD RE SPECIFlED ON THE 

C DATA CARP. THE INITIAL POPULATION VECTOR WILL RE COMPtITED FROM THE 

c STABLE AGE I,ISTRIRlITION. 

C FOLLOWED BY THE CPRRELATIONS OF FECUNDITY RATES, BETWEEN FECUNDITY 

C AMI, SURVIVAL RATES, AND AMONG SURVIVAL RATES. 

C.....THE COVARIANCE MATRIX (IW SINGLE SURSCRIPT FORM) IS FIRST REA@ IWTO 

C WORK VECTOR W, THEN REPACKED IN STARDARD UPPER TRIANGULAR FORM IN 

C.....COVAR. 

C.....LOAP 

c.....IF ~0~2x7, THE COVAR DATA CARD SHOULD CONTAIN THE ?*N VARIANCES, 

DIMENSION X(1) ,BX(l) ,PX(l) ,COVAR(l) ,W(1) 

10 FORMAT (8Fln.6) 

11 FORMAT (8F10.3) 

12 FCIRMAT (11) 

13 FORMAT (I?) 

14 FORMAT (15) 

READ (1,13) N 

REAP (1,13) PJT 

READ (1,14) NR 

READ (1,12) TOP1 

READ (1,13) TOP? 

READ (1,1?) LUN 

READ (1,14) ISEEP 

READ (1,13) KK 

READ (1,lO) RLl 

READ (1,lO) I4IDTH2 

REAn (1,lO) RL 

REAP (1,in) WIDTHI 

REAP (i,in) RL~ 
READ (1,in) RU 

IF (Inpi) irm,iono,iozn inon READ (1,ii) (x(I),I=i,ri) 

GO TO 1041) 

in40 NZ=N*Z 

REAP (1,lO) (BX(I),I=l,N) 

REAP (1,lO) (PX(I),I=l,N) 

1020 REAP (1,ll) X(l) 

NCV=N* (N2+1) 

C.. ... COVARIANCE VECTOR 

IF (IOPZ) inm,~~m,in~n inm READ (1,in) (w(r),~=i,~cv) 

in80 REAP (1,in) (w(I),I=~,N~) 

GO TO 3110 

READ (1,lO) C1 

READ (1,lO) C7 

IF (ABS(C1)-1.) 40,4O,90 

READ (1,in) c3 

fin IF (ABs(cz)-i.) 5n,50,90 

Fin IF (ms(c3)-1.) ~in,m,9n 

60 no 70 I=l,N? 

IF (M(1)) 8O,70,70 

GO TO 100 

70 CONTINUE 

80 WRITE (LiiN,m) 

GO TO 3znn 

90 WRITE (LuN,w) 

 XIHA HAVE ABSOLUTE VALUE GREATER THAN 1.0~ ) 

GO TO 3znn 

89 FORMAT (/I ERROR -- VARIAWCES IN COVARIANCE MATRIX CANNOT'/ 

f9X'HAVE NEGATIVE VALUES ') 

99 FORMAT (/I ERROR -- ELEMENTS QF CORRELATION MATRIX CANNOT'/ 

C.....REPACK 

100 MID=NCV-N*(N+1)/3 

NDX=NCV+l 

K=NCV 

I=N2 

M=l 

DO 12031) NNN=l,NCV 

NDXzNDX- 1 

IF (NDX-K) iinn,ii6o,iinn 

iino IF (NDX-MID) iiio,iiin,iizo 

iiin IF (VDX-L) ii5n,ii40,ii4n 

iim W(WOX>=CZ*SCRT(W( I )*W(J) 

iiw w(NDx)=c~*s~RT(w(I)*w(J)) 

1120 W( NDX)=C3*SORT( W( I)*W( J)) 

GO TO izon 

GO TO 1200 

GO TO i2nn 

1161) W(NDX)=W( I) 

L=K-N 

M=M+1 

K=K-M 

I=I-1 
51 

J=N3+1 

1200 3=J-1 

C.....ZERO-MEAN CASES 

2mn DO 3100 I=~,N 

IF ( BX ( I - . onooooi) 2030,703(1,7010 

znin IF (px(I)-.oaooooi) mn,zozn,2ion 

mn IR=I+N 

GO TO 3040 

204n DO 2070 J=I,FP 

IF (J-IR) ~050,2050,2060 

2030 IR=I 

2050 I J= (J -1) *N2- ( (5-3)*5) /2-J+I R 

GO TO 2070 

2060 IJ=( IR-l)*N2-( (IR-3)*IR)/Z-IR+J 

z07n w ( I J)=o. 0 

2100 CONTINUE 

ziin NDX=O 

IF (IR-I-N) ZO~O,ZIOO,~O 

C.. . . .REPACK AGAIN TO STANDARD llPPER TRIANGULAR FORM 

00 2140 I=l,N? 

DO 2120 J=I,N? 

NDX=NDX+l 

I ,I= ( J *J - J ) /2+1 

COVAR( IJ)=W(NDX) 

2120 CONTINUE 

2140 CONTINUE 

END 

znn RETURN 

SllRROlJT I NE EREAL (RNX , ENX, N, NT, ER, ERT) 

C*****COMPUTES THE REALIZED FACTOR OF INCREASE FOR AN INITIAL POPULATION 

c VECTOR RNX ANI, A FINAL VECTOR ENX OF ORPER N OVER NT TIME STEPS. 

DIMENSION RNX(1) ,ENX(l) 

TI=O.O 

TF=O. 0 

TI=T I+RNX ( I ) 

TF=TF+ENX( I) 

ERT=TF/T I 

ER=ERT**( 1. O/FLOAT( NT) ) 

RETURN 

END 

DO loo0 I=~,P 

ioon CONTINIJE 

S\lRROUTINE RPI N (RX ,PX,W ,N) 

C*****LOADS VECTORS BX AND PX INTO W. 

RIMENSION RX(1) ,PX(l),W(l) 

NDX=N 

NDX=NDX+l 

DO inoo I=I,N 

W( I)=RX( I) 

W(NDX)=PX( I) 

inno CONTINUE 

RETURV 
52 

E NI, 

SURROUTINE RPOUT (W,BT,PT,N) 

DIMENSION W(1) ,RT(l) ,PT(1) 

NDX=N 

ND X= N D X+ 1 

RT(I)=W( I) 

PT( I)=W(NDX) 

1QQO CONTINUE 

RETURN 

ENT, 

C*****\INLOAPS THE CONTENTS OF W INTO BT AFJD PT. 

no inm I=I,M 
RECENT TECHNICAL MEMORANDUMS 

Copies of this and other NOAA Technical Memorandums are available from the National Technical 

Information Service, 5285 Port Royal Road, Springfield, VA 22167. Paper copies vary in price. 

Microfiche copies cost $3 50. Recent issues of NOAA Technical Memorandums from the NMFS 

Southwest Fisheries Center are listed below 

NOAA TM NMFS SLYFC 18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

Environmental data contouring program EDMAPZ 

L EBER 

(April 1982) 

The relationship between changes in gross reproductive rate and the 

current rate of increase for some simple aqe structured models 

T POLACHECK 

(May 1982) 

Testing methods of estimating range and bearing to cetaceans 

T D SMITH 

(1982) 

"An annotated bibliography of the ecology of co-occurring tunas 

(Katsuwonus pelamis. Thunnus albacares) and dolphins (Stenella 

attenuata. Stenella longirostris and Delphinus delphis in the 

eastern tropical Pacific" 

aboard the R/V D S Jordan 

S D HAWES 

(November 1982) 

Structured flotsam as fish aggregating devices 

R S SHOMURAandW M MATSUMOTO 

(November 1982) 

Abundance estimation of dolphin stocks involved in the eastern 

tropical Pacific yellowfin tuna fishery determined from aerial and 

ship surveys to 1979. 

R. S. HOLT and J. E. POWERS 

(November 1982) 

Revised update and retrieval system for the CalCOFl oceanographic 

L. EBER and N. WlLEY 

(December 1982) 

A preliminary study of dolphin release procedures using model purse 

D. B. HOLTS and J. M. COE 

(December 1982) 

"Possible effects of sampling biases on reproduction rate estimates 

TOM POLACHECK 

(January 1983) 

"Report of porpoise experiment testing detection of on-track schools 

RENNIE S. HOLT 

(February 1983) 

data file. 

seines. 

for porpoise in the eastern tropical Pacific." 

(pet dots), March 7-April 5, 1981 "