model
	{

		for (i in 1 : N) {

			for(t in 1:Nyr){

				TFR[i,t]       ~ dnorm(m[i,t], lambda.f[i,t])
				m[i,t]         <- Fpost[i] + pF[i,t] *(Fpre[i] - Fpost[i]) 
				logit(pF[i,t]) <- beta[i]*(year[t] - tau[i]) 
				lambda.f[i,t]  <- women[i,t] * prec.f

			}

			Fpre[i]  <- a1 + d1 * educ70[i] + e1[i]
			Fpost[i] <- a2 + d2 * educ00[i] + e2[i]
			t10[i]   <- a3 + d3 * educ70[i] + e3[i]
			t90[i]   <- a4 + d4 * educ70[i] + e4[i]

                        beta[i]  <- -2*log(9)/(t90[i]-t10[i])
                        tau[i]   <- (t10[i] + t90[i])/2
	                    
                        start[i]    <- t10[i]
                        duration[i] <- t90[i]-t10[i]

                        # Forecasting the fertility rate in year 2005  

                        logit(pF2005[i]) <-  beta[i]*(4.5 - tau[i])    
                        m2005[i]         <-  Fpost[i] + pF2005[i] *(Fpre[i] - Fpost[i])  
                        aux2005[i]       <- (m2005[i]*(women[i,5]/100))/tot

                        # Forecasting the fertility rate in year 2010  

                        logit(pF2010[i]) <-  beta[i]*(5.0 - tau[i])    
                        m2010[i]         <-  Fpost[i] + pF2010[i] *(Fpre[i] - Fpost[i])  
                        aux2010[i]       <- (m2010[i]*(women[i,5]/100))/tot
                }
	
		# total number of women in 2000 dividied by 100
		# To use as the weight denominator in the forecasting
		tot <- sum(women[ , 5])/100
		
		# Forecasting TFR in Brazil in 2010
		tfr2010 <- sum(aux2010[ ])
		tfr2005 <- sum(aux2005[ ])
			
		# CAR prior distribution for spatial random effects: 
		e1[1:N] ~ car.normal(adj[], weights[], num[], lambda1)
		e2[1:N] ~ car.normal(adj[], weights[], num[], lambda2)
		e3[1:N] ~ car.normal(adj[], weights[], num[], lambda3)
		e4[1:N] ~ car.normal(adj[], weights[], num[], lambda4)
	
		for(k in 1:sumNumNeigh) {
			weights[k] <- 1
	}
	

#-- PRIORS

	prec.f    ~ dgamma(0.01, 0.01)

	a1 ~ dflat() 
	a2 ~ dflat() 
	a3 ~ dflat() 
	a4 ~ dflat() 
		
	d1 ~ dnorm(0.0, 0.001)
	d2 ~ dnorm(0.0, 0.001)
	d3 ~ dnorm(0.0, 0.001)
	d4 ~ dnorm(0.0, 0.001)
		
	lambda1 ~ dgamma(0.01, 0.01) 
	lambda2 ~ dgamma(0.01, 0.01) 
	lambda3 ~ dgamma(0.01, 0.01)
	lambda4 ~ dgamma(0.01, 0.01)		

		
        sigma1 <- 1/sqrt(lambda1)
        sigma2 <- 1/sqrt(lambda2)
        sigma3 <- 1/sqrt(lambda3)
        sigma4 <- 1/sqrt(lambda4)

        sigma.f <- 1/sqrt(prec.f)
			
	}