##############################################################################################################################
##############################################################################################################################
##R CODE FOR ANNUAL POPULATION ESTIMATES BY AGE AND SEX, GIVEN LITTLE DATA (SINGLE COUNTY)
##
##EDDIE HUNSINGER, JANUARY 2011 (LAST UPDATED SEPTEMBER 5, 2011)
##http://www.demog.berkeley.edu/~eddieh/
##edynivn@gmail.com
##
##A WRITEUP WITH FULL EXPLANATION IS FORTHCOMING, BUT HERE IS A LINK TO THE INITIAL IDEA I HAD (EDITS SINCE): http://www.google.com/reader/item/tag:google.com,2005:reader/item/50ae883a25bee913
##THE ESTIMATES USE/REQUIRE THE COUNTY COUNTS BY AGE AND SEX FROM THE LAST TWO CENSUSES, A GENERIC GROSS MIGRATION PROFILE, A GENERIC LIFE TABLE AND ANNUAL COUNTS FOR BIRTHS AND NET MIGRATION
##
##TO RUN THE CODE, YOU SHOULD BE ABLE TO SIMPLY SELECT ALL AND PASTE INTO R
##IDEALLY, THE USER WOULD CONTROL (WITH ITERATIVE PROPORTIONAL FITTING (http://www.demog.berkeley.edu/~eddieh/datafitting.html)) COUNTY ESTIMATES MADE THIS WAY UP TO MORE CAREFULLY MADE STATE ESTIMATES
##ALSO, THE USER WOULD MAKE MODIFICATIONS OR DISAGGREGATIONS FOR ANY PERCEIVED PROBLEMS, SUCH AS ESTIMATION OF PARTICULAR SPECIAL POPULATIONS
##THIS CODE/PROCESS MIGHT BE PARTICULARLY USEFUL FOR VERY SMALL AREAS, PERHAPS CONTROLLED UP TO COUNTIES
##
##IT IS CURRENTLY SET TO RUN Douglas COUNTY FEMALES, BUT YOU CAN CHANGE/REPLACE (JUST REPLACE ALL) THE AREA NAME TO ANY COLORADO COUNTY, AND/OR TO MALES (M)
##CODE TO RUN ALL OF THE COUNTIES AT ONCE IS AVAILABLE AT http://www.demog.berkeley.edu/~eddieh/AppliedDemographyToolbox/AgeSexCountyEstimates/AgeSexCountyCode_All.txt
##
##IF YOU WOULD LIKE TO USE, SHARE OR REPRODUCE ANY INFORMATION OR IDEAS FROM THIS WORK, BE SURE TO CITE THE SOURCE
##THERE IS NO WARRANTY FOR THIS CODE
##
##PLEASE NOTE THAT I'VE FOUND PROBLEMS WITH THIS CODE (SEE STEP 11), AND IT IS A WORK IN PROGRESS
##############################################################################################################################
##############################################################################################################################

##############################################################################################################################
##############################################################################################################################
##Some more notes: (1) In this proposed way, migration is applied to age groups, and not cohorts--this may seem unrealistic, but if, alternatively, 
##migration were applied to the cohorts (x-(10-n)), I think areas with special populations (such as college campuses) would be very poorly accounted, 
##and, also, I'm hoping that controlling to more carefully done state estimates will smooth the problems of applying migration to the age groups. 
##(2) This work may be nothing new at all-- I believe it is the same as a "plus-minus adjustment", and it parts that are like "Hamilton-Perry", and many 
##other residual estimates works. 
##############################################################################################################################
##############################################################################################################################

##DIMENSIONS
SIZE<-91

##############################################################################################################################
#STEP 1: Doing one sex at a time, read in census county populations by age
##############################################################################################################################

K<-read.table(file="http://www.demog.berkeley.edu/~eddieh/AppliedDemographyToolbox/AgeSexCountyEstimates/AgeSexCounty_90and00.csv",header=TRUE,sep=",")
KCensus1<-K$Douglas_F_1990
KCensus2<-K$Douglas_F_2000

##############################################################################################################################
#STEP 2: Subtract any "special population" (population by age that you want to hold fixed (such as college students) or manage separately 
#from the 1990 and 2000 county populations.
##############################################################################################################################

Adjust<-read.table(file="http://www.demog.berkeley.edu/~eddieh/AppliedDemographyToolbox/AgeSexCountyEstimates/SpecialPopAdjustments_90and00.csv",header=TRUE,sep=",")
Adjust1<-Adjust$Douglas_F_1990
Adjust2<-Adjust$Douglas_F_2000
KBase1<-KCensus1-Adjust1
KBase2<-KCensus2-Adjust2

##############################################################################################################################
#STEP 3: Read in data on total annual births and set the fraction female at birth and fraction male at birth
##############################################################################################################################

Births<-read.table(file="http://www.demog.berkeley.edu/~eddieh/AppliedDemographyToolbox/AgeSexCountyEstimates/Births_90sand00s.csv",header=TRUE,sep=",")
Births<-Births$Douglas
ffab<-.4886
fmab<-1-ffab
fab<-ffab

##############################################################################################################################
#STEP 4: Read in a generic survival profile (for instance, I just use one for the US, 2000) and calculate Sx for the projection matrix
##############################################################################################################################

Mortality<-read.table(file="http://www.demog.berkeley.edu/~eddieh/AppliedDemographyToolbox/AgeSexCountyEstimates/lx_00.csv",header=TRUE,sep=",")
lx<-Mortality$lx2000
Lx<-array(0,c(SIZE))
for (i in 1:SIZE){Lx[i]<-.5*(lx[i]+lx[i+1])}
Sx<-array(0,SIZE-1)
for (i in 1:SIZE-1){Sx[i]<-(Lx[i+1]/Lx[i])}
Survival<-array(0,c(SIZE,SIZE))
Survival[,]<- rbind(0,cbind(diag(Sx),0))

##############################################################################################################################
#STEP 5: Read in a generic gross migration profile (sums to one) (for instance, I just use one for Colorado, 1995-2000)
##############################################################################################################################

GrossMigration<-read.table(file="http://www.demog.berkeley.edu/~eddieh/AppliedDemographyToolbox/AgeSexCountyEstimates/GrossMigrationProfile_95to00.csv",header=TRUE,sep=",")
Gx<-GrossMigration$GrossMigProf

##############################################################################################################################
#STEP 6: Read in data on total annual net migration
##############################################################################################################################

NetMigration<-read.table(file="http://www.demog.berkeley.edu/~eddieh/AppliedDemographyToolbox/AgeSexCountyEstimates/NetMigration_90sand00s.csv",header=TRUE,sep=",")
NetMigration<-NetMigration$Douglas

##############################################################################################################################
#STEP 7: Make arrays to hold the data for each year
##############################################################################################################################

K20<-K19<-K18<-K17<-K16<-K15<-K14<-K13<-K12<-K11<-K10<-K9<-K8<-K7<-K6<-K5<-K4<-K3<-K2<-K1<-K0<-array(0,c(SIZE,1))

##############################################################################################################################
#STEP 8: Simply age the population one year forward (including rates of survival), add births and repeat
##############################################################################################################################

Births0<-array(0,c(SIZE,1))
Births0[1]<-Births[1]
K1<-(Survival%*%KBase1)+(Births0)*fab

Births1<-array(0,c(SIZE,1))
Births1[1]<-Births[2]
K2<-(Survival%*%K1)+(Births1)*fab

Births2<-array(0,c(SIZE,1))
Births2[1]<-Births[3]
K3<-(Survival%*%K2)+(Births2)*fab

Births3<-array(0,c(SIZE,1))
Births3[1]<-Births[4]
K4<-(Survival%*%K3)+(Births3)*fab

Births4<-array(0,c(SIZE,1))
Births4[1]<-Births[5]
K5<-(Survival%*%K4)+(Births4)*fab

Births5<-array(0,c(SIZE,1))
Births5[1]<-Births[6]
K6<-(Survival%*%K5)+(Births5)*fab

Births6<-array(0,c(SIZE,1))
Births6[1]<-Births[7]
K7<-(Survival%*%K6)+(Births6)*fab

Births7<-array(0,c(SIZE,1))
Births7[1]<-Births[8]
K8<-(Survival%*%K7)+(Births7)*fab

Births8<-array(0,c(SIZE,1))
Births8[1]<-Births[9]
K9<-(Survival%*%K8)+(Births8)*fab

Births9<-array(0,c(SIZE,1))
Births9[1]<-Births[10]
K10<-(Survival%*%K9)+(Births9)*fab

##############################################################################################################################
#STEP 9: Calculate net migration by age (simply the second census population minus the 10th projection step)
##############################################################################################################################

Nx<-KBase2-K10

##############################################################################################################################
#STEP 10: Standardize the net migration by proportionally adjusting it with the generic gross migration profile
##############################################################################################################################

NxStandard<-Nx-Gx*sum(Nx)

##############################################################################################################################
#STEP 11: Apply the net migration profile to each intercensal year by adding the sum of net migration through the specified intercensal year 
#***multiplied by a proportion of net migration (naive .5 is used on this example, but is not recommended for all cases) for the sub-population (males)***, 
#multiplied by Gx plus Nxstandard times n/10 (the standard profile is interpolated) to Kx[t+n]temp 
#
#Note in a previous version I did the following to get the sub-population ratio:
#***multiplied by the sum of the residually-calculated net migration (Nx) divided by the sum of the annual net migration data you read in (this ratio provides 
#a simple proportional adjustment by sex)***
#I found this creates significant problems for a lot of cases and determined that I either need to integrate the sub-groups or use something empirical.
##############################################################################################################################

K1<-K1+(sum(NetMigration[1])*.5)*Gx+NxStandard*.1
K2<-K2+(sum(NetMigration[1:2])*.5)*Gx+NxStandard*.2
K3<-K3+(sum(NetMigration[1:3])*.5)*Gx+NxStandard*.3
K4<-K4+(sum(NetMigration[1:4])*.5)*Gx+NxStandard*.4
K5<-K5+(sum(NetMigration[1:5])*.5)*Gx+NxStandard*.5
K6<-K6+(sum(NetMigration[1:6])*.5)*Gx+NxStandard*.6
K7<-K7+(sum(NetMigration[1:7])*.5)*Gx+NxStandard*.7
K8<-K8+(sum(NetMigration[1:8])*.5)*Gx+NxStandard*.8
K9<-K9+(sum(NetMigration[1:9])*.5)*Gx+NxStandard*.9
K10<-K10+(sum(NetMigration[1:10])*.5)*Gx+NxStandard*1

##############################################################################################################################
#STEP 12: Repeat STEP 8 for the second decade (post-censal)
##############################################################################################################################

Births10<-array(0,c(SIZE,1))
Births10[1]<-Births[11]
K11<-(Survival%*%KBase2)+(Births10)*fab

Births11<-array(0,c(SIZE,1))
Births11[1]<-Births[12]
K12<-(Survival%*%K11)+(Births11)*fab

Births12<-array(0,c(SIZE,1))
Births12[1]<-Births[13]
K13<-(Survival%*%K12)+(Births12)*fab

Births13<-array(0,c(SIZE,1))
Births13[1]<-Births[14]
K14<-(Survival%*%K13)+(Births13)*fab

Births14<-array(0,c(SIZE,1))
Births14[1]<-Births[15]
K15<-(Survival%*%K14)+(Births14)*fab

Births15<-array(0,c(SIZE,1))
Births15[1]<-Births[16]
K16<-(Survival%*%K15)+(Births15)*fab

Births16<-array(0,c(SIZE,1))
Births16[1]<-Births[17]
K17<-(Survival%*%K16)+(Births16)*fab

Births17<-array(0,c(SIZE,1))
Births17[1]<-Births[18]
K18<-(Survival%*%K17)+(Births17)*fab

Births18<-array(0,c(SIZE,1))
Births18[1]<-Births[19]
K19<-(Survival%*%K18)+(Births18)*fab

Births19<-array(0,c(SIZE,1))
Births19[1]<-Births[20]
K20<-(Survival%*%K19)+(Births19)*fab

##############################################################################################################################
#STEP 13: Repeat STEP 11 for the second decade (post-censal)
##############################################################################################################################

K11<-K11+(sum(NetMigration[11])*.5)*Gx+NxStandard*.1
K12<-K12+(sum(NetMigration[11:12])*.5)*Gx+NxStandard*.2
K13<-K13+(sum(NetMigration[11:13])*.5)*Gx+NxStandard*.3
K14<-K14+(sum(NetMigration[11:14])*.5)*Gx+NxStandard*.4
K15<-K15+(sum(NetMigration[11:15])*.5)*Gx+NxStandard*.5
K16<-K16+(sum(NetMigration[11:16])*.5)*Gx+NxStandard*.6
K17<-K17+(sum(NetMigration[11:17])*.5)*Gx+NxStandard*.7
K18<-K18+(sum(NetMigration[11:18])*.5)*Gx+NxStandard*.8
K19<-K19+(sum(NetMigration[11:19])*.5)*Gx+NxStandard*.9
K20<-K20+(sum(NetMigration[11:20])*.5)*Gx+NxStandard*1

##############################################################################################################################
#STEP 14: Replace any negative numbers with zeroes (if/when you adjust the estimates up to other estimates, be sure to make the 
#zeroes "flexible" (eg make them .001 so that they can proportionally adjust)
##############################################################################################################################

K10[K10<0]<-K9[K9<0]<-K8[K8<0]<-K7[K7<0]<-K6[K6<0]<-K5[K5<0]<-K4[K4<0]<-K3[K3<0]<-K2[K2<0]<-K1[K1<0]<-0
K20[K20<0]<-K19[K19<0]<-K18[K18<0]<-K17[K17<0]<-K16[K16<0]<-K15[K15<0]<-K14[K14<0]<-K13[K13<0]<-K12[K12<0]<-K11[K11<0]<-0

##############################################################################################################################
#STEP 15: Do what you like with the data formatting and read it out to a csv
##############################################################################################################################

K0<-KCensus1
K1<-K1+Adjust1
K2<-K2+Adjust1
K3<-K3+Adjust1
K4<-K4+Adjust1
K5<-K5+Adjust1
K6<-K6+Adjust1
K7<-K7+Adjust1
K8<-K8+Adjust1
K9<-K9+Adjust1
K10<-KCensus2
K11<-K11+Adjust2
K12<-K12+Adjust2
K13<-K13+Adjust2
K14<-K14+Adjust2
K15<-K15+Adjust2
K16<-K16+Adjust2
K17<-K17+Adjust2
K18<-K18+Adjust2
K19<-K19+Adjust2
K20<-K20+Adjust2

plot(K20,type="l",col="red")
points(KCensus2,type="l",lwd=5)
points(K10,type="l",col="red",lwd=5)
points(K9,type="l",col="blue")
points(K8,type="l",col="green")
points(K7,type="l",col="purple")
points(K6,type="l",col="red")
points(K5,type="l",col="blue")
points(K4,type="l",col="green")
points(K3,type="l",col="purple")
points(K2,type="l",col="red")
points(K1,type="l",col="blue")
points(KCensus1,type="l",lwd=5)
points(K20,type="l",col="blue",lwd=5)
points(K19,type="l",col="red")
points(K18,type="l",col="green")
points(K17,type="l",col="purple")
points(K16,type="l",col="red")
points(K15,type="l",col="blue")
points(K14,type="l",col="green")
points(K13,type="l",col="purple")
points(K12,type="l",col="red")
points(K11,type="l",col="blue")

#write.table(K20, file="G:/###/###.csv", sep=",")