b01_Demographics; Demographics.bas
Attribute VB_Name = "b01_Demographics"
Option Explicit
Dim num As Double
Dim n As Long, i As Long, a_count As Long
Dim ptemp() As Double ' For alignment
Dim a_pool() As Long, f_pool() As Long
'!*****************************************************
'! This module handles demographics
'!*****************************************************
Public Sub new_population()
'! Main module demographics
Printdok "new_population"
Dim i As Long, h As Long
' Updating some variables
Call UpdateVariables
' Mortality
Call Mortality
'Adoption
Call adopt
' Migration
Call migration
' Fertility
Call Fertility
' Young people leaving home
Call leave_home
' Bilda par
Call cohab
' Split pair
Call Split_Pair
' Early retirement
' Note: if running SESIM in "BabyBoom mode" another model is used to simulate
' intry into disability pension
If get_scalefactor("BabyBoom_Active") <> 1 Then
Call DisabilityPension
Else
Call early_retire
End If
' Rehab (exiting early retirement status)
Call rehab
' Do checks of consistence
Call do_checks("after demographics")
' Update descriptive demographics
Call demograf_stat
' Write life events access file (if enabled in options)
lifehist.write_now
Printdok "-- new_population ready"
End Sub
'*****************************************************************************
' Subroutine for prediction of death for each individual. The results can be
' subjected to alignment and/or variance reduction. Alignment makes sure that
' either the average expected hazard within each alignment stratum equals the
' prespecified level or that the expected number of deaths within each alignment
' stratum equals the prespecified level.
' Variance reduction techniques reduces the Monte Carlo variance introduced
' by the simulation so that the simulated deathcounts equals, or is very
' close to, the expected count or the alignment count.
'*****************************************************************************
Private Sub Mortality()
Printdok "Mortality" '*** print to log-file
Call stock_sub 'calculations for death-risk, used in pension model
Dim alignment As Integer '*** 1 - alignment on, 0 - alignment off
Dim var_reduction As Integer '*** 1 - variance reduction on, 0 - v.red. off
Dim model As Integer '*** 1 - full model, 2 - simple model (Statistics Sweden)
Dim alignment_type As Integer '*** 1 - by rates, 2 - by numbers
Dim pi_death() As Double '*** death-risks
Dim wanted(1 To 20, 1 To 2) As Double '*** wanted number of deaths (age-group / sex)
Dim expected(1 To 20, 1 To 2) As Double '*** expected number of deaths (age-group / sex)
Dim run_time_mortality As Double
Dim counts(1 To 20, 1 To 2) As Long, scaleparms(1 To 20, 1 To 2) As Double
Dim agegrp As Integer, sex As Integer, randomnr As Long
Dim scaleparms2(1 To 20, 1 To 2) As Double, counts2(1 To 20, 1 To 2) As Long
randomnr = Int(Rnd * 100000)
' read run-time scale factor for mortality (default value: 1)
run_time_mortality = get_scalefactor("mortality")
'*** Parameters that control alignment, variance reduction and choice of
'*** model. The parameters will later be controlled properly through the interface
alignment = get_scalefactor("mort_align")
var_reduction = get_scalefactor("mort_varred")
model = get_scalefactor("mort_mod")
alignment_type = get_scalefactor("alignment_type")
'*** default values
If alignment < 0 Or alignment > 1 Then alignment = 1
If var_reduction < 0 Or var_reduction > 1 Then var_reduction = 1
If model < 1 Or model > 2 Then model = 1
If alignment_type < 1 Or alignment_type > 2 Then alignment_type = 1
ReDim pi_death(1 To m_icount) As Double
'*** Calculate risks for each individual
Call calc_deathprobs(model, pi_death)
'*** Calculate excess mortality for individuals with functional disorder
'*** (only when the Baby Boom project is active)
' NOTE: Due to the increased level of mortality caused by the rescaling of
' the hazards have to be aligned to the hazards (or counts) of Statistics Sweden
If get_scalefactor("BabyBoom_Active") <> 1 Then
Call ExcessMortality(pi_death)
alignment = 1
var_reduction = 1
End If
If var_reduction = 0 Then
If alignment = 0 Then
'*** No alignment or variance reduction. The calculated probabilities are
'*** applied directly through Monte-Carlo simulation.
'*** scale parameters equals the run time mortality (default = 1)
For agegrp = 1 To 20
For sex = 1 To 2
scaleparms(agegrp, sex) = run_time_mortality
Next sex
Next agegrp
Call rand_mort(1, scaleparms, pi_death, counts)
Else '*** alignment = 1
'*** NOTE: Alignment without variance reduction is here achieved by simple
'*** proportional scaling of the probabilities. If this case should be more
'*** frequently used in the future the routine should be rewritten to make
'*** the adjustment on another scale, e.g. the logit scale. Using proportional
'*** adjustment large errors can occur for probabilities close to zero or one.
'*** TP 021016
'*** Calculate wanted counts for alignment and/or variance reduction
Call calc_wanted_mortality(1, wanted, pi_death, alignment_type)
'*** Calculate expected counts for alignment and/or variance reduction
Call calc_wanted_mortality(0, expected, pi_death, 1)
'*** Calculate scale parameters using ratio between wanted and
'*** expected counts (agegroup and sex)
For agegrp = 1 To 20
For sex = 1 To 2
' update wanted count with run time scaling
wanted(agegrp, sex) = wanted(agegrp, sex) * run_time_mortality
If expected(agegrp, sex) > 0 Then
scaleparms(agegrp, sex) = CDbl(wanted(agegrp, sex) / expected(agegrp, sex))
Else
scaleparms(agegrp, sex) = 1
End If
Next sex
Next agegrp
Call rand_mort(1, scaleparms, pi_death, counts)
End If '*** alignment = 1
Else '*** var_reduction = 1
'*** Variance reduction with or without alignment. Without alignment in this
'*** case rather means alignment to the expected number of events given the
'*** used probabilities and individuals at risk. With alignment means alignment
'*** to exogenously given counts (or rates).
'*** Variance reduction is achieved using the "raising the tripwire"-algorithm
'*** (see documentation). TP 021016
Dim i As Long, j As Long
Dim idx As Double, totcount As Double, rand() As Double
Dim totwant As Double, counter(1 To 20, 1 To 2) As Double, rand2 As Double
Dim diff1 As Double, diff2 As Double, pi As Double
Dim sortvec() As Single
Dim age As Integer
ReDim sortvec(1 To m_icount, 1 To 2), rand(1 To m_icount)
'*** Calculate wanted counts for alignment and/or variance reduction
Call calc_wanted_mortality(alignment, wanted, pi_death, alignment_type)
'*** Draw uniform random numbers
Call RANUNI(m_icount, rand(1), model_time + base_year + random * Rnd)
'*** Create matrix containing the logit-transformed variate and individual number
For i = 1 To m_icount
rand2 = mini(0.9999999, maxi(0.0000001, rand(i)))
pi = mini(0.9999999, maxi(0.0000001, pi_death(i)))
sortvec(i, 1) = logit(rand2) - logit(pi)
sortvec(i, 2) = i_indnr(i)
Next
'*** Free memory
Erase pi_death, rand
'*** Sort matrix by logit variate
'*** NOTE: Too large matrices causes errors in the call to SORTMATRIX due
'*** to stack overflows.
If m_icount <= 100000 Then
Call SORTMATRIX(m_icount, 2, sortvec(1, 1))
Else
Call quicksort(sortvec, 1)
End If
'*** Loop across the sorted matrix and select individuals until
'*** each alignment stratum is filled
'*** NOTE2: the method can give a larger number of events than what is
'*** wanted if the absolute deviation between the number of calculated
'*** and wanted events is minimized.
totcount = 0
For i = 1 To m_icount
'*** individual index number
idx = indnr2index(sortvec(i, 2))
agegrp = mini(Int(i_age(idx) / 5) + 1, 20)
mark_i(idx) = 0
'*** Decide the optimal stopping point - before or after one
'*** additional event?
diff1 = Abs(counter(agegrp, i_sex(idx)) - wanted(agegrp, i_sex(idx)))
diff2 = Abs(counter(agegrp, i_sex(idx)) - wanted(agegrp, i_sex(idx)) + m_weight)
'*** Check if alignment stratum is "full"
If counter(agegrp, i_sex(idx)) < wanted(agegrp, i_sex(idx)) Then
'*** Is the alignment stratum filled after the next event?
If ((counter(agegrp, i_sex(idx)) + m_weight) < wanted(agegrp, i_sex(idx))) Then
mark_i(idx) = 1
'*** NOTE: individuals living abroad is not counted for alignment
If i_abroad(idx) = 0 Then
counter(agegrp, i_sex(idx)) = counter(agegrp, i_sex(idx)) + m_weight
totcount = totcount + m_weight
End If
Else
'*** Adding one event gives occurred > wanted but is allowed if
'*** abs(occurred - wanted) decreases
If diff1 > diff2 Then
mark_i(idx) = 1
'*** NOTE: individuals living abroad is not counted for alignment
If i_abroad(idx) = 0 Then
counter(agegrp, i_sex(idx)) = counter(agegrp, i_sex(idx)) + m_weight
totcount = totcount + m_weight
End If
End If
End If
End If
Next
'*** Calculate the total number of wanted deaths
For i = 1 To 20
For j = 1 To 2
totwant = totwant + wanted(i, j)
Next
Next
m_dead = totcount
status "Dead: " & m_dead & " (" & totwant & ")"
'*** Optional debug output - set runtime parameter
If get_scalefactor("debug_mortality") <> 1 Then
Dim filenr As Integer
Dim debugfile As String
filenr = FreeFile
debugfile = sesimpath & "\tempdata\debug_mort.txt"
If model_time = 1 And Dir(debugfile) <> "" Then
Open debugfile For Output As #filenr
Else
Open debugfile For Append As #filenr
End If
If model_time = 1 Then Print #filenr, _
"year, dead, age, sex, edlevel, pgi, civstat, status, ftp_64"
For i = 1 To m_icount
Write #filenr, base_year + model_time, mark_i(i), i_age(i), _
i_sex(i), i_edlevel(i), i_pgi(i), i_civ_stat(i), _
i_status(i), i_ftp_64(i)
Next
Close #filenr
End If
Call delete_individuals
End If '*** var_reduction = 1
Call dödsrisk_sub 'calculations for death-risk, used in pension model
End Sub
'*****************************************************************************
' Subroutine for calculation of hazards. The current version supports two
' mortality models as specified below.
' Arguments:
' model(IN): 1 - one-year hazards (by sex and age) projected by Statistics
' Sweden, excess mortality for the early retired lowers the general
' mortality so that the total mortality (for living in Sweden) is
' unchanged.
' 2 - model estimated on LINDA 1997
' probs(OUT): array containing calculated hazards
'*****************************************************************************
Private Sub calc_deathprobs(ByVal model As Integer, ByRef probs() As Double)
Dim year As Integer, age As Integer, count_dead As Long
' For each age: 1 - not early retired, men; 2 - not early retired, women
' 3 - early retired, men, 4 - early retired, women
Dim count(17 To 64, 1 To 4) As Long
' Hazards for all, for early retired and for non early-retired
Dim haz As Double, haz_ftp As Double, haz_non_ftp As Double
Dim n_all As Long, n_ftp As Long, n_not_ftp As Long
If model = 2 Then
'*** Risks according to projection from Statistics Sweden with adjustment
'*** for excess mortality of the early retired. NOTE - this leaves the
'*** total mortality rate (by age and sex) unaffected.
year = base_year + model_time
' Get population counts for not early retired and early retired respectively
For i = 1 To m_icount
'*** NOTE: only individuals living in Sweden are counted for alignment
If i_age(i) >= 17 And i_age(i) <= 64 And i_abroad(i) = 0 Then
If i_status(i) <> 4 Then
count(i_age(i), i_sex(i)) = count(i_age(i), i_sex(i)) + 1
Else
count(i_age(i), 2 + i_sex(i)) = count(i_age(i), 2 + i_sex(i)) + 1
End If
End If
Next
ReDim mark_i(1 To UBound(mark_i)) As Long '*** erase previous markers
For i = 1 To m_icount
age = mini(i_age(i), 105)
haz = parm_death(mini(2110, year), age, i_sex(i))
probs(i) = haz
' Scale down probabilites for not early retired
If age >= 17 And age <= 64 Then
haz_ftp = parm_death_er(age, i_sex(i))
n_all = count(i_age(i), i_sex(i)) + count(i_age(i), 2 + i_sex(i))
n_ftp = count(i_age(i), 2 + i_sex(i))
n_not_ftp = count(i_age(i), i_sex(i))
If n_not_ftp > 0 Then
haz_non_ftp = (n_all * haz - n_ftp * haz_ftp) / n_not_ftp
Else
haz_non_ftp = haz_ftp
End If
If i_status(i) = 4 Then
probs(i) = haz_ftp
Else
probs(i) = haz_non_ftp
End If
End If '*** age >= 17 And age <= 64
Next i
Else
'********************************************************************
'*** Risks as estimated from LINDA - see documentation at
'*** N:\FI_E4\Sesim\Estimering\Mortalitet
'********************************************************************
Dim civstat As Integer, ftp As Integer, bkon As Integer
Dim pgi_q() As Double
Dim xbeta As Double, pgi As Double, ed As Double
'*** calculate quintile one for PGI (30 - 64 years) and total pensions (65 -)
Dim pgitemp() As Double
ReDim pgitemp(1 To m_icount) As Double
Dim counter1 As Long
counter1 = 0
' prepare temporary vectors
For i = 1 To m_icount
Select Case (i_age(i))
Case 30 To 64
counter1 = counter1 + 1
pgitemp(counter1) = i_pgi(i)
Case Else
End Select
Next i
' remove redundant elements
ReDim Preserve pgitemp(1 To counter1) As Double
pgi_q = arr_Percentile(pgitemp, 20, 40, 60, 80)
year = base_year + model_time
For i = 1 To m_icount
Select Case (i_age(i))
'*** For individuals aged 0 to 29, hazards from the demographic projection
'*** of Statistics Sweden (1999) are used
Case Is <= 29
age = mini(i_age(i), 105)
probs(i) = parm_death(mini(2110, year), age, i_sex(i))
Case 30 To 64
'*** code variables
If i_civ_stat(i) = 0 Then civstat = 1 Else civstat = -1
If i_status(i) = 4 Then ftp = -1 Else ftp = 1
If i_sex(i) = 1 Then bkon = 1 Else bkon = -1
If i_pgi(i) <= pgi_q(1, 2) Then pgi = 0.395777289
If pgi_q(1, 2) < i_pgi(i) And i_pgi(i) <= pgi_q(2, 2) Then pgi = 0.25577989
If pgi_q(2, 2) < i_pgi(i) And i_pgi(i) <= pgi_q(3, 2) Then pgi = -0.033832204
If pgi_q(3, 2) < i_pgi(i) And i_pgi(i) <= pgi_q(4, 2) Then pgi = -0.22024363
If i_pgi(i) > pgi_q(4, 2) Then pgi = -0.397481345
xbeta = -9.087658517 + _
bkon * 0.37214669 + _
i_age(i) * 0.070821524 + _
ftp * -1.451930375 + _
i_age(i) * ftp * 0.016452455 + _
pgi + _
civstat * 0.237121176
probs(i) = 1 / (1 + Exp(-xbeta))
'*** Risks from SCB are used for emigrants since they have no income
If i_abroad(i) = 1 Then
age = mini(i_age(i), 105)
probs(i) = parm_death(mini(2110, year), age, i_sex(i))
End If
Case Is >= 65
'*** code variables
If i_civ_stat(i) = 0 Then civstat = 1 Else civstat = -1
If i_ftp_64(i) = 1 Then ftp = -1 Else ftp = 1
If i_sex(i) = 1 Then bkon = 1 Else bkon = -1
If i_edlevel(i) = 0 Then ed = 0.149320967
If i_edlevel(i) = 1 Then ed = -0.006917504
If i_edlevel(i) = 2 Then ed = -0.142403462
xbeta = -11.83489131 + _
i_age(i) * 0.112117111 + _
ftp * -1.420891233 + _
bkon * 1.210972038 + _
civstat * 0.968908827 + _
i_age(i) * civstat * -0.010775078 + _
i_age(i) * ftp * 0.015966928 + _
i_age(i) * bkon * -0.011559863 + _
ed
probs(i) = 1 / (1 + Exp(-xbeta))
'*** Risks from SCB are used for emigrants since they have no income
If i_abroad(i) = 1 Then
age = mini(i_age(i), 105)
probs(i) = parm_death(mini(2110, year), age, i_sex(i))
End If
End Select
Next i
End If
End Sub
'*****************************************************************************
' Sub calc_wanted_mortality calculates the wanted number of deaths. If
' alignment is used the wanted number of deaths returned is either the expected
' number of deaths with respect to the risks used or the number of deaths
' as specified by the forecast of Statistics Sweden depending on the value of
' the parameter atype.
' If variance reduction is used without alignment implicit alignment towards
' the expected count is made. The expected counts are calculated using the
' current population within each alignment stratum and the age-specific hazards
' given by Statistics Sweden.
' Arguments:
' align(IN): 0 - no alignment, 1 - alignment
' wanted(OUT): matrix by agegroup and sex containing the calculated
' number of wanted counts
' probs(IN): array containing calculated hazards
' atype(IN): 1 - alignment by rates, 2 - alignment by numbers
'*****************************************************************************
Private Sub calc_wanted_mortality(ByVal align As Integer, ByRef wanted() As Double, _
ByRef probs() As Double, ByVal atype As Integer)
Dim agegroup As Integer, sex As Integer
If align = 0 Then '*** No alignment used: wanted = expected ***
Dim i As Long
For i = 1 To m_icount
agegroup = mini(Int(i_age(i) / 5) + 1, 20) ' 0-4, 5-9, ..., 95-
'*** NOTE: only individuals living in Sweden are counted for alignment
If i_abroad(i) = 0 Then
wanted(agegroup, i_sex(i)) = wanted(agegroup, i_sex(i)) + probs(i) * m_weight
End If
Next
Else '*** Alignment used: wanted = external data ***
Dim year As Integer, age As Integer
year = base_year + model_time
'*** Alignment by numbers
If atype = 2 Then
'*** Get the current year from parm_align_mortality
For age = 1 To 20
For sex = 1 To 2
wanted(age, sex) = parm_align_mortality(year, age, sex)
Next
Next
'*** Alignment by rates
Else
For i = 1 To m_icount
agegroup = mini(Int(i_age(i) / 5) + 1, 20) ' 0-4, 5-9, ..., 95-
age = mini(i_age(i), 105)
'*** NOTE: only individuals living in Sweden are counted for alignment
If i_abroad(i) = 0 Then
wanted(agegroup, i_sex(i)) = wanted(agegroup, i_sex(i)) + _
parm_death(mini(2110, year), age, i_sex(i)) * m_weight
End If
Next
End If
End If
'*** Round counts to integer values
For agegroup = 1 To 20
For sex = 1 To 2
wanted(agegroup, sex) = round(wanted(agegroup, sex), 0)
Next
Next
End Sub
'*****************************************************************************
' Function for randomization of deaths due to calculated probabilities.
' Arguments:
' active(IN): 0 - no action taken, 1 - executes deaths
' scalefactor(IN): matrix containing scalefactors used for
' alignment/variance reduction
' probs(IN): array containing calculated hazards
' counts(OUT): matrix containing counts of deaths
'*****************************************************************************
Private Sub rand_mort(ByVal active As Integer, scaleparms() As Double, _
ByRef probs() As Double, ByRef counts() As Long)
Printdok "die" ' print to log-file
Dim i As Long, flag As Integer, count_dead As Long, count_dead_total As Long
Dim agegrp As Integer, sex As Integer, tempcount(1 To 20, 1 To 2) As Double
count_dead = 0
count_dead_total = 0
flag = 0
For i = 1 To m_icount
agegrp = mini(Int(i_age(i) / 5) + 1, 20)
mark_i(i) = 0
If Rnd < probs(i) * scaleparms(agegrp, i_sex(i)) Then
mark_i(i) = 1
flag = 1
'*** total count of death including emigrants
count_dead_total = count_dead_total + 1
'*** NOTE: only individuals living in Sweden are counted for alignment
If i_abroad(i) = 0 Then
tempcount(agegrp, i_sex(i)) = tempcount(agegrp, i_sex(i)) + _
m_weight '*** count by group
'*** macro count of deaths
'*** NOTE: only persons living in Sweden
count_dead = count_dead + 1
End If
End If
Next i
'*** Round macro group-counts to integer values
For agegrp = 1 To 20
For sex = 1 To 2
counts(agegrp, sex) = round(tempcount(agegrp, sex), 0)
Next sex
Next agegrp
If active = 1 Then status "Dead: " & count_dead * m_weight & " (" & _
count_dead_total * m_weight & ")"
m_dead = count_dead * m_weight
If active = 1 And flag = 1 Then
If controlcenter.chk2SaveAccessdb.value = 1 Then
status "Saving dead peoples last record..."
Call MDIForm1.menu_writeaccess_Click
End If
Call död_sub 'calculations for death-risk, used in pension model
Call delete_individuals 'delete individuals from data vectors
End If
End Sub
'**************************************************************************
'*** Sub adopt() distributes orphans to new households. Note that orphans
'*** older than 17 years are automatically assigned the headship of a new
'*** household.
'**************************************************************************
Private Sub adopt()
Printdok "Adopt"
Dim i As Long, n As Long, f_count As Long, min_age As Long, max_age As Long
Dim SAMBO As Integer, emam As Integer, epap As Integer, age As Long
Dim age2 As Long, m As Long, f As Long, limit As Long, diff As Long
ReDim a_pool(m_hcount), f_pool(m_hcount), ptemp(m_hcount)
Dim indnr As Long
Dim flag As Integer
a_count = 0
f_count = 0
Call clear_marks
flag = 0
' Select orphan housholds.
' Children 18 years of age or older are not adopted but forms new households.
For i = 1 To m_hcount
If h_n_child(hhnr2index(h_hhnr(i))) > 0 And h_n_adults(hhnr2index(h_hhnr(i))) = 0 Then
' Mark all children that are not to be adopted.
min_age = 18
indnr = h_first_indnr(hhnr2index(h_hhnr(i)))
Do Until indnr = 0
If i_age(indnr2index(indnr)) < min_age Then min_age = i_age(indnr2index(indnr))
If i_age(indnr2index(indnr)) >= 18 Then
mark_i(indnr2index(indnr)) = 1
flag = 1
End If
indnr = i_next_indnr(indnr2index(indnr))
Loop
' Are there any children that are to be adopted?
If min_age < 18 Then
a_count = a_count + 1
a_pool(a_count) = h_hhnr(i)
End If
End If
Next
' Create new households for the old orphans using the same routine as for
' youngsters leaving home
If flag = 1 Then Call leave_home_doit
'*** Adoptions are handled only if there are orphans to adopt...
If a_count > 0 Then
'Select possible families and assign "probabilities"
SAMBO = 0
emam = 0
epap = 0
f_count = 0
For i = 1 To m_hcount
If h_n_adults(hhnr2index(h_hhnr(i))) > 0 And h_abroad(i) = 0 Then
Call get_malefemale_indnr(h_hhnr(i), m, f)
'Age criteria
If m = 0 Then
max_age = i_age(indnr2index(f))
min_age = i_age(indnr2index(f))
ElseIf f = 0 Then
max_age = i_age(indnr2index(m))
min_age = i_age(indnr2index(m))
Else
max_age = maxi(i_age(indnr2index(m)), i_age(indnr2index(f)))
min_age = mini(i_age(indnr2index(m)), i_age(indnr2index(f)))
End If
If min_age = 0 Then min_age = max_age
If max_age < 51 And min_age > 24 Then
f_count = f_count + 1
mark_h(hhnr2index(h_hhnr(i))) = 1
'Family type and age
'Age refers to the female (male if no adult female in hh)
If h_n_adults(i) = 2 Then 'Cohabitants
SAMBO = 1
age = i_age(indnr2index(f))
Else
If f = 0 Then 'Single male
epap = 1
age = i_age(indnr2index(m))
Else 'Single female
emam = 1
age = i_age(indnr2index(f))
End If
End If
age2 = age * age
'Predicted - actual number of children
diff = -6.90615 _
- 0.846521 * emam _
- 1.328339 * epap _
+ 0.472703 * age _
- 0.00643 * age2 _
- h_n_child(i) _
+ gauss(0, 0.91058)
'Store differences in vector ptemp
ptemp(f_count) = diff
End If '*** max_age < 51 And min_age > 24
End If
Next
limit = findlimit(ptemp(), a_count, 0)
f_count = 0
n = 0
' Create pool of adopting households (by household identity numbers)
For i = 1 To m_hcount
If mark_h(i) = 1 Then f_count = f_count + 1
' Compare difference to limit
If mark_h(i) = 1 And ptemp(f_count) < limit Then mark_h(i) = 0 ' Not selected
' A happy family who is going to adopt a child
If mark_h(i) = 1 And h_abroad(i) = 0 Then
n = n + 1
f_pool(n) = h_hhnr(i)
End If
Next
'Nu vet vi vilka barn som behöver föräldrar och vilka familjer som skall adoptera, Matcha!!
status "Adopt " & a_count * m_weight 'Ger antal familjer, inte antal barn
Call match_families
End If '*** a_count > 0
End Sub
Private Sub match_families()
Dim f As Long, indnr As Long, firstind As Long
Call clear_marks
For f = 1 To a_count
'lifehist.write_hist a_pool(f), "Adopted by " & f_pool(f)
' Mark orphan hh for delete
mark_h(hhnr2index(a_pool(f))) = 1
'lifehist.write_hist a_pool(f), "Delete hh (adopted)"
' new hhnr to all members in childs hh from new family
firstind = h_first_indnr(hhnr2index(a_pool(f)))
indnr = firstind
Do Until indnr = 0
i_hhnr(indnr2index(indnr)) = f_pool(f)
lifehist.write_hist indnr, "New hhnr (adopted)"
indnr = i_next_indnr(indnr2index(indnr))
Loop
' update hh-info
h_n_child(hhnr2index(f_pool(f))) = _
h_n_child(hhnr2index(f_pool(f))) + h_n_child(hhnr2index(a_pool(f)))
h_size(hhnr2index(f_pool(f))) = _
h_size(hhnr2index(f_pool(f))) + h_n_child(hhnr2index(a_pool(f)))
' set last ind in family hh poiting to first ind in child hh
i_next_indnr(indnr2index(get_hhlast_indnr(f_pool(f)))) = firstind
Next
' Resten av koden tagen direkt från match_couples
' Time to delete households
Dim j As Long, i As Long
j = m_hcount
For i = 1 To m_hcount
If mark_h(i) = 1 Then
hhnr2index(h_hhnr(i)) = 0
j = j - 1
End If
Next
Call shift_up_marked_hh
m_hcount = j
For i = 1 To m_hcount
hhnr2index(h_hhnr(i)) = i
' Clear vector
mark_h(i) = 0
Next
' Keeping the household vectors packed to prevent errors
Call dyn_vect_h(m_hcount)
Erase a_pool, f_pool, ptemp
End Sub
'*****************************************************************************
'*** Sub Fertility is used to randomize newborn children to mothers.
'*** Alignment and variance reduction can be used in the same fashion as in
'*** sub mortality.
'*** For fertility the alignment total is disaggregated into five-year
'*** groups with respect to womens ages.
'*****************************************************************************
Private Sub Fertility()
Printdok "Fertility" '*** print to log-file
Dim alignment As Integer '*** 1 - alignment on, 0 - alignment off
Dim var_reduction As Integer '*** 1 - variance reduction on, 0 - v.red. off
Dim pi_fertility() As Double '*** death-risks
Dim wanted As Double '*** wanted number of deaths
Dim expected As Double '*** expected number of deaths
Dim run_time_fertility As Double
Dim count As Long, scaleparm As Double
Dim agegrp As Integer, sex As Integer, randomnr As Long
Dim scaleparm2 As Double, count2 As Long
randomnr = Int(Rnd * 100000)
' read run-time scale factor for mortality
run_time_fertility = get_scalefactor("fertility")
'*** Parameters that control alignment, variance reduction and choice of
'*** model. The parameters will be controlled through the user interface later...
alignment = get_scalefactor("fert_align")
var_reduction = get_scalefactor("fert_varred")
' default values
If alignment < 0 Or alignment > 1 Then alignment = 1
If var_reduction < 0 Or var_reduction > 1 Then var_reduction = 1
ReDim pi_fertility(1 To m_icount) As Double
'*** Calculate risks for each individual
Call calc_fertilityprobs(pi_fertility)
If var_reduction = 0 Then
If alignment = 0 Then
scaleparm = run_time_fertility
Call rand_fertility(1, scaleparm, pi_fertility, count)
Else '*** alignment = 1
'*** Calculate wanted counts for alignment and/or variance reduction
Call calc_wanted_fertility(1, wanted, pi_fertility)
'*** Calculate expected counts for alignment and/or variance reduction
Call calc_wanted_fertility(0, expected, pi_fertility)
'*** Calculate scale parameter using ratio between wanted and
'*** expected count
' wanted count updated with run time scaling factor
wanted = wanted * run_time_fertility
If expected > 0 Then
scaleparm = CDbl(wanted / expected)
Else
scaleparm = 1
End If
Call rand_fertility(1, scaleparm, pi_fertility, count)
End If '*** alignment = 1
Else '*** var_reduction = 1
'*** Variance reduction with or without alignment. Without alignment in this
'*** case rather means alignment to the expected number of events given the
'*** used probabilities and individuals at risk. With alignment means alignment
'*** to exogenously given counts (or rates).
'*** Variance reduction is achieved using the "raising the tripwire"-algorithm
'*** (see documentation). TP 021016
'*** Age groups for alignment: 15-19, 20-24,..., 45-49
Dim i As Long, j As Long
Dim idx As Double, totcount As Double, rand() As Double
Dim totwant As Double, counter As Double, rand2 As Double
Dim diff1 As Double, diff2 As Double, pi As Double
Dim sortvec() As Single
Dim age As Integer
Dim tfrcalc(14 To 50, 1 To 2) As Long
Dim ageidx As Long
ReDim sortvec(1 To m_icount, 1 To 2), rand(1 To m_icount)
'*** Reset the matrix used for calculation of TFR
For i = 14 To 50
tfrcalc(i, 1) = 0
tfrcalc(i, 2) = 0
Next
'*** Calculate wanted counts for alignment and/or variance reduction
Call calc_wanted_fertility(alignment, wanted, pi_fertility)
'*** Draw uniform random numbers
Call RANUNI(m_icount, rand(1), model_time + base_year + random * Rnd)
'*** Create matrix containing the logit-transformed variate and individual number
For i = 1 To m_icount
rand2 = mini(0.9999999, maxi(0.0000001, rand(i)))
pi = mini(0.9999999, maxi(0.0000001, pi_fertility(i)))
sortvec(i, 1) = logit(rand2) - logit(pi)
sortvec(i, 2) = i_indnr(i)
Next
'*** Free memory
Erase pi_fertility, rand
'*** Sort matrix by logit variate
'*** NOTE: Too large matrices causes errors in the call to SORTMATRIX due
'*** to stack overflows.
If m_icount <= 100000 Then
Call SORTMATRIX(m_icount, 2, sortvec(1, 1))
Else
Call quicksort(sortvec, 1)
End If
'*** Loop across the sorted matrix and select individuals until
'*** each alignment stratum is filled
'*** NOTE2: the method can give a larger number of events than what is
'*** wanted if the absolute deviation between the number of calculated
'*** and wanted events is minimized.
counter = 0
For i = 1 To m_icount
'*** individual index number
idx = indnr2index(sortvec(i, 2))
mark_i(idx) = 0
'*** Check if individual is at risk
If i_sex(idx) = 2 And i_age(idx) >= 18 And i_age(idx) <= 49 And _
i_bvux(idx) And i_abroad(idx) = 0 Then
'*** Calculate individuals at risk
tfrcalc(i_age(idx), 1) = tfrcalc(i_age(idx), 1) + 1
'*** Decide the optimal stopping point - before or after one
'*** additional event?
diff1 = Abs(counter - wanted)
diff2 = Abs(counter - wanted + m_weight)
'*** Check if alignment stratum is "full"
If counter < wanted Then
'*** Is the alignment stratum filled after the next event?
If (counter + m_weight) < wanted Then
mark_i(idx) = 1
counter = counter + m_weight
'*** Accumulate number of born
tfrcalc(i_age(idx), 2) = tfrcalc(i_age(idx), 2) + 1
Else
'*** Adding one event gives occurred > wanted but is allowed if
'*** abs(occurred - wanted) decreases
If diff1 > diff2 Then
mark_i(idx) = 1
counter = counter + m_weight
'*** Accumulate number of born
tfrcalc(i_age(idx), 2) = tfrcalc(i_age(idx), 2) + 1
End If
End If
End If '*** Stratum full?
End If '*** Individual at risk?
Next
'*** Calculate TFR
m_TFR = 0
For i = 14 To 50
If tfrcalc(i, 1) > 0 Then
m_TFR = m_TFR + CDbl(tfrcalc(i, 2)) / CDbl(tfrcalc(i, 1))
End If
Next
m_born = counter
status "Born: " & m_born & " (" & wanted & ")"
Call add_born
End If '*** var_reduction = 1
End Sub
'*****************************************************************************
'*** calc_fertilityprobs calculates probabilities of giving birth during the
'*** year based on estimated logistic regression models.
'*** For further documentation see N:\FI_E4\Sesim\Estimering\Fertility
'***
'*** Arguments:
'*** probs(OUT): calculated probabilities
'*****************************************************************************
Private Sub calc_fertilityprobs(ByRef probs() As Double)
Dim i As Long, indnr As Long, count As Long
Dim n_child As Integer, agechild As Integer
Dim temp_pgi() As Double, pgi_q() As Double
Dim da1 As Integer, da2 As Integer, da3 As Integer, da4 As Integer, da5 As Integer
Dim da6 As Integer, da7 As Integer, q1 As Integer, q2 As Integer, q3 As Integer
Dim ed1 As Integer, ed2 As Integer, bforv As Integer
Dim xbeta As Double
' calculate quartiles for PGI (used as covariate)
' NOTE: only for individuals at risk
count = 0
ReDim temp_pgi(1 To m_icount) As Double
For i = 1 To m_icount
If i_sex(i) = 2 And i_age(i) >= 18 And i_age(i) <= 49 And _
i_bvux(i) = 1 And i_abroad(i) = 0 Then
count = count + 1
temp_pgi(count) = i_pgi(i)
End If
Next i
ReDim Preserve temp_pgi(1 To count) As Double 'remove empty space
pgi_q = arr_Percentile(temp_pgi, 25, 50, 75)
For i = 1 To m_icount
' Probabilities are calculated for "adult" women in ages 18-49 living in sweden
If i_sex(i) = 2 And i_age(i) >= 18 And i_age(i) <= 49 And _
i_bvux(i) = 1 And i_abroad(i) = 0 Then
' determine the number of children in ages 0-15 and the age of the
' youngest child
agechild = 15
n_child = 0
If h_size(hhnr2index(i_hhnr(i))) > h_n_adults(hhnr2index(i_hhnr(i))) Then
indnr = h_first_indnr(hhnr2index(i_hhnr(i)))
Do Until indnr = 0
If i_bvux(indnr2index(indnr)) = 0 Then
n_child = n_child + 1
If i_age(indnr2index(indnr)) < agechild Then _
agechild = i_age(indnr2index(indnr))
End If
indnr = i_next_indnr(indnr2index(indnr))
Loop
End If
'*** general covariates
da1 = CInt(((i_age(i) >= 18) And (i_age(i) <= 19))) * -1
da2 = CInt(((i_age(i) >= 20) And (i_age(i) <= 24))) * -1
da3 = CInt(((i_age(i) >= 25) And (i_age(i) <= 29))) * -1
da4 = CInt(((i_age(i) >= 30) And (i_age(i) <= 34))) * -1
da5 = CInt(((i_age(i) >= 35) And (i_age(i) <= 39))) * -1
da6 = CInt(((i_age(i) >= 40) And (i_age(i) <= 44))) * -1
da7 = CInt(((i_age(i) >= 45) And (i_age(i) <= 49))) * -1
q1 = CInt(i_pgi(i) <= pgi_q(1, 2)) * -1
q2 = CInt(i_pgi(i) > pgi_q(1, 2) And i_pgi(i) <= pgi_q(2, 2)) * -1
q3 = CInt(i_pgi(i) > pgi_q(2, 2) And i_pgi(i) <= pgi_q(3, 2)) * -1
ed1 = CInt(i_edlevel(i) = 0) * -1
ed2 = CInt(i_edlevel(i) = 1) * -1
bforv = CInt(i_status1(i) = 8 And i_inc_taxable1(i) > 0) * -1
' temporary solution until i_status1 and i_income1 are created in base data
If model_time = 1 Then
If i_status(i) = 8 Then bforv = 1 Else bforv = 0
End If
'*** Different models for different parities
Select Case n_child
' first child
Case 0
xbeta = -17.4674 + _
i_age(i) * 0.8516 + _
CInt(i_age(i)) * CInt(i_age(i)) * -0.016 + _
i_civ_stat(i) * 4.3097 + _
q1 * 1.0481 + _
q2 * 1.396 + _
q3 * 0.5528 + _
bforv * 0.5792 + _
ed1 * 0.2798 + _
ed2 * -0.0174
' second child
Case 1
xbeta = -5.810965 + _
da2 * 0.356915 + _
da3 * 0.702668 + _
da4 * 1.148072 + _
da5 * 1.331858 + _
da6 * 0.720743 + _
i_civ_stat(i) * 2.322846 + _
q1 * 0.991257 + _
q2 * 0.963913 + _
q3 * 0.203133 + _
bforv * 0.472624 + _
ed1 * -0.524592 + _
ed2 * -0.335218 + _
agechild * 1.236023 + _
agechild * i_age(i) * -0.037041
' third child
Case 2
xbeta = -5.7038 + _
da2 * 1.017409 + _
da3 * 0.67553 + _
da4 * 0.772821 + _
da5 * 1.020751 + _
da6 * 0.840791 + _
i_civ_stat(i) * 0.958956 + _
q1 * 1.206739 + _
q2 * 0.857448 + _
q3 * 0.008772 + _
bforv * 0.312307 + _
ed1 * -0.293166 + _
ed2 * -0.488417 + _
agechild * 1.110564 + _
i_age(i) * agechild * -0.032289
' fourth child
Case 3
xbeta = -16.650261 + _
da2 * 11.820241 + _
da3 * 12.029253 + _
da4 * 12.068464 + _
da5 * 12.243513 + _
da6 * 12.56333 + _
i_civ_stat(i) * 0.249647 + _
q1 * 0.923402 + _
q2 * 0.678916 + _
q3 * 0.092358 + _
bforv * -0.062157 + _
ed1 * 0.463188 + _
ed2 * 0.19574 + _
agechild * 1.148298 + _
i_age(i) * agechild * -0.03442
Case Else
xbeta = -500
End Select
'*** Alignment of fertility: the model seems to seriously underpredict the level of
'*** fertility when not aligned. Due to this a simple alignment constant is used to
'*** boost the levels temporarily until the fertility model is reestimated. The
'*** fertility model is aligned using forecasts from Statistics Sweden until 2050
'*** and the below alignment adjusts for the discrepancy in fertility levels between
'*** 2051 and the following years.
'*** TP060516
xbeta = xbeta + 0.8
probs(i) = 1 / (1 + Exp(-xbeta))
' NOTE: a maximum of four births per woman is allowed!!
If n_child >= 4 Then probs(i) = 0
Else
' If not "adult" women in ages 18-49 living in sweden
probs(i) = 0
End If
Next i
End Sub
'*****************************************************************************
'*** calc_wanted_fertility calculates the wanted / expected number of births
'*** during the year.
'***
'*** Arguments:
'*** align(IN): 0 - wanted count equals expected count, 1 - wanted count as
'*** projected by Statistics Sweden
'*** wanted(OUT): calculated counts
'*** probs(IN): probabilities for calculation of expected counts
'*****************************************************************************
Private Sub calc_wanted_fertility(ByVal align As Integer, ByRef wanted As Double, _
ByRef probs() As Double)
Dim sex As Integer
'*** No alignment used: wanted = expected ***
If align = 0 Then
Dim i As Long
For i = 1 To m_icount
'*** NOTE: only "adult" women living in Sweden in ages 18 - 49 are counted for alignment
If i_sex(i) = 2 And i_abroad(i) = 0 And i_age(i) >= 18 And i_age(i) <= 49 _
And i_bvux(i) = 1 Then
wanted = wanted + probs(i) * m_weight
End If
Next i
'*** Round count to integer values
wanted = round(wanted, 0)
'*** Alignment used: wanted = external data ***
Else
Dim year As Integer
year = base_year + model_time
wanted = parm_align_fertility(mini(2110, year))
'*** To prevent a drop in the TFR the level is raised from 2051 and
'*** onwards ???? ThP 070309, läser in serier tom 2110
'If base_year + model_time >= 2051 Then
' wanted = wanted * 1.03
'End If
End If
End Sub
'*****************************************************************************
'*** rand_fertility randomizes births to specific women based their estimated
'*** probabilities.
'***
'*** Arguments:
'*** active(IN): 0 - calculates the number of randomized individuals, no action
'*** taken, 1 - randomization performed
'*** scaleparms(IN): factors for proportional scaling of probabilities (by agegroup)
'*** probs(IN): estimated probabilities
'*** counts(OUT): calculated counts (by womens agegroup) of newborn
'*****************************************************************************
Private Sub rand_fertility(ByVal active As Integer, scaleparm As Double, _
ByRef probs() As Double, ByRef count As Long)
Printdok "fertility" ' print to log-file
Dim i As Long, flag As Integer, count_born As Long
Dim tempcount As Double
Dim tfrcalc(14 To 50, 1 To 2) As Long
count_born = 0
flag = 0
For i = 1 To m_icount
mark_i(i) = 0
' Only "adult" women in ages 18-49 living in sweden are at risk
If i_abroad(i) = 0 And i_age(i) >= 18 And i_age(i) <= 49 And i_sex(i) = 2 And _
i_bvux(i) = 1 Then
'*** Calculate # of individuals at risk
tfrcalc(i_age(i), 1) = tfrcalc(i_age(i), 1) + 1
If Rnd < probs(i) * scaleparm Then
mark_i(i) = 1
flag = 1
count_born = count_born + 1 '*** macro count
tempcount = tempcount + m_weight
'*** Accumulate # of individuals born per mothers age
tfrcalc(i_age(i), 2) = tfrcalc(i_age(i), 2) + 1
End If
End If
Next i
'*** Round macro count to integer values
count = round(tempcount, 0)
If active = 1 Then status "Born: " & count_born * m_weight
m_born = count_born * m_weight
'*** Calculate tfr
m_TFR = 0
For i = 14 To 50
If tfrcalc(i, 1) > 0 Then
m_TFR = m_TFR + CDbl(tfrcalc(i, 2)) / CDbl(tfrcalc(i, 1))
End If
Next
If active = 1 And flag = 1 Then Call add_born
End Sub
'*********************************************************************************
' Module for decision to leave home and form new household.
' The probabilities are estimated from HEK 99 for males and females separately in
' age 18-29.
' Alignment to the expected value is implemented in order to reduce the Monte-Carlo
' variation.
'
' NOTE: These estimates refer to the marginal probability of an individual being
' adult (having left his parents home) at a certain age given the covariates. The
' simulations are however based on the hazards, i.e. the probability of moving out
' at age t given that the individual has not moved out until time t-1. Using the
' marginal probailities in place of the conditional probabilities (hazards) the
' model will overestimate the number of individuals leaving home.
' Below the calculated marginal probabilites are adjusted using ad-hoc factors
' to give the correct marginal distribution of adults by age and sex.
' The marginal probabilities could be recalculated to the proper conditional
' probabilities using P_t_cond = (P_t - P_t-1) / (1 - P_t-1), where P_t denotes
' the estimated marginal probability at age t. This is based on the assumption that
' the other covariates does not change between years.
'*********************************************************************************
Private Sub leave_home()
'! Youngsters leaving home
Printdok "leave_home"
Dim i As Long, n_atrisk As Long
Dim flag As Integer
Dim haz() As Double, expval As Double, scalefactor As Double, atrisk_indnr() As Double
ReDim atrisk_indnr(1 To m_icount), haz(1 To m_icount)
'Male parameters
Const m0 As Single = -28.01799 'Intercept
Const m1 As Single = 2.16731 'bald
Const m2 As Single = -0.03854 'bald*bald
Const m3 As Single = -1.69635 'Low education
Const m4 As Single = -1.83159 'Medium education
Const m5 As Single = 0.05349 'bald*Low education
Const m6 As Single = 0.06191 'bald*Medium education
Const m7 As Single = -0.20749 '1 if study allowance > 0, else 0
Const m8 As Single = 0.25834 '1 if taxable income > 0, else 0
Const m9 As Single = -0.23805 'Nationality: 1 if Swedish
'Female parameters
Const f0 As Single = -23.07892 'Intercept
Const f1 As Single = 1.71356 'bald
Const f2 As Single = -0.02869 'bald*bald
Const f3 As Single = -2.07506 'Low education
Const f4 As Single = -1.56884 'Medium education
Const f5 As Single = 0.1033 'bald*Low education
Const f6 As Single = 0.06794 'bald*Medium education
Const f7 As Single = -0.21878 '1 if study allowance > 0, else 0
Const f8 As Single = 1.00927 '1 if taxable income > 0, else 0
Const f9 As Single = -0.3211 'Nationality: 1 if Swedish
Dim v1, v2, v3, v4, v5, v6, v7, v8, v9, v10 As Double
Dim ex As Double
' Calculate expected value and create a pool of individuals at risk
expval = 0
n_atrisk = 0
For i = 1 To m_icount
If i_bvux(i) = 0 And h_n_adults(hhnr2index(i_hhnr(i))) > 0 And i_age(i) >= 18 Then
'*** At ages 30 and above probabilities are calculated using age 29
'*** TP021010
v1 = mini(29, i_age(i))
v2 = v1 * v1
If i_edlevel(i) = 0 Then v3 = 1 Else v3 = 0
If i_edlevel(i) = 1 Then v4 = 1 Else v4 = 0
v5 = v1 * v3
v6 = v1 * v4
If i_trf_study(i) > 0 Then v7 = 1 Else v7 = 0
If i_inc_taxable1(i) > 0 Then v8 = 1 Else v8 = 0
If i_born_abroad(i) = 0 Then v9 = 0 Else v9 = 1
If i_sex(i) - 1 = 0 Then
ex = Exp(m0 + m1 * v1 + m2 * v2 + m3 * v3 + m4 * v4 + m5 * v5 + m6 * v6 + _
m7 * v7 + m8 * v8 + m9 * v9)
Else
ex = Exp(f0 + f1 * v1 + f2 * v2 + f3 * v3 + f4 * v4 + f5 * v5 + f6 * v6 + _
f7 * v7 + f8 * v8 + f9 * v9)
End If
n_atrisk = n_atrisk + 1
atrisk_indnr(n_atrisk) = i_indnr(i)
'*** Adjustment of probabilities TP021010
If i_sex(i) = 1 Then
haz(n_atrisk) = (1 / (1 + Exp(-ex))) * 0.3
Else
haz(n_atrisk) = (1 / (1 + Exp(-ex))) * 0.27
End If
' Print_to_file "debug_leave_home.txt", "N", i_age(i), i_sex(i), haz(n_atrisk)
expval = expval + haz(n_atrisk) * m_weight
End If
Next
'*** Delete redundant elements in vector haz
ReDim Preserve haz(1 To n_atrisk)
Dim sortvek() As Double, randvek() As Double
ReDim sortvek(1 To n_atrisk), randvek(1 To n_atrisk)
' Draw uniform random numbers
Call RANUNI(n_atrisk, randvek(1), model_time + base_year + random * Rnd)
' Calculate the logit-transformed variable for sorting-based alignment
Dim idx As Long
For i = 1 To n_atrisk
idx = indnr2index(atrisk_indnr(i))
' NOTE: findlimit returns threshold for the K LARGEST values. Hence, we
' calculate the negative of the standard sorting variable
sortvek(i) = -(logit(randvek(i)) - logit(haz(i)))
Next
' Find the limit for returning the K largest values
' NOTE: naive simulation would use limit = 0!!
Dim limit As Double
limit = findlimit(sortvek, CLng(expval / m_weight), 0)
Dim counter As Long
flag = 0
counter = 0
Call clear_marks
' Mark the chosen individuals
For i = 1 To n_atrisk
If sortvek(i) > limit Then
counter = counter + 1
mark_i(indnr2index(atrisk_indnr(i))) = 1
mark_h(hhnr2index(i_hhnr(indnr2index(atrisk_indnr(i))))) = 1
flag = 1
End If
' End the loop - we're done
If counter = CLng(expval / m_weight) Then Exit For
Next
status "Leaving home: " & counter * m_weight
If flag = 1 Then Call leave_home_doit
End Sub
'*************************************************************************************
'*** COMMENT?
'*************************************************************************************
Public Sub migration()
Printdok "migration"
Call emigration
Call immigration
m_netmigration = m_immigrated - m_emigrated
End Sub
Public Sub immigration()
'The following routine has the following purpose:
' 1) Assume gross immigration, X, read from Excel-sheet
' 2) Simulate number of re-immigrating Swedes, Xs. This is a fixed number
' which by definition equals the number of emigrating Swedes, Ys.
' 3) Clone X-Xs=Xi non-swedish immigrantes from previously immigrated non-swedish households
'--------------------------------------------------------------------------
'Next steg checks whether the household resides abroad
'If so is the case a random number is generated and compared to the hazard
'of returning to Sweden conditional on number of years abroad
Dim n_wantedS As Long
Dim n_wantedI As Long
Dim n_actualS As Long
Dim n_actualI As Long
Dim year As Integer
Dim hh As Long
Dim hhi As Long
Dim inr As Long
Dim index As Long
Dim rand As Double
Dim pred As Double
Dim tid As Integer
year = model_time + base_year
Dim scale_par As Double
scale_par = 1 / m_weight
n_wantedS = scale_par * n_im(mini(2110, year), 1) 'Migrationsdata börjar på 1998, är det ett problem???
n_wantedI = scale_par * n_im(mini(2110, year), 2) 'Migrationsdata börjar på 1998, är det ett problem???
n_actualS = 0 'Antal personer som återimmigrerat i snurran nedan.
n_actualI = 0 'Antal personer som återimmigrerat i snurran nedan.
Dim hazard(0 To 100) As Double
Dim survival(0 To 100) As Double
Dim x As Integer
Const intercept As Double = 2.3144
Const alfa As Double = 1.5
Const epsilon As Double = 5
hazard(0) = 0
hazard(1) = Exp(-(intercept - 0.27707))
hazard(2) = Exp(-(intercept - 1.48956))
hazard(3) = Exp(-(intercept - 1.45464))
hazard(4) = Exp(-(intercept - 1.29084))
hazard(5) = Exp(-(intercept - 1.01211))
hazard(6) = Exp(-(intercept - 0.79006))
hazard(7) = Exp(-(intercept - 0.39685))
'Using a pareto distribution to assign hazard values where time>7
For x = 8 To 100
' hazard(x) = hazard(0)
hazard(x) = alfa * (epsilon ^ alfa) * x ^ (-alfa - 1)
' Debug.Print hazard(x)
Next
'Calculating survival curve
survival(0) = 1
For i = 1 To 100
survival(i) = survival(i - 1) * (1 - hazard(i))
Next
ReDim dist_n_wantedS(1 To 2, 1 To 100) As Double 'Number of returning Swedes wanted distributed on time abroad
Dim sum_n_wantedS As Double
Dim last_hazard As Integer
Dim rest_sum As Integer
sum_n_wantedS = 0
last_hazard = 0
rest_sum = 0
'Calculating the pdf from the hazard/survival information and the number
'of individuals to return after time=i conditional on n_wantedS. Periods
'where the number of individuals is rounded to 0 are clustered in rest_sum.
For i = 1 To 100
dist_n_wantedS(1, i) = hazard(i) * survival(i - 1) * n_wantedS
sum_n_wantedS = sum_n_wantedS + dist_n_wantedS(1, i)
If dist_n_wantedS(1, i) < 0.5 Then
last_hazard = i
rest_sum = n_wantedS - sum_n_wantedS
Exit For
End If
Next
ReDim Preserve dist_n_wantedS(1 To 2, 1 To last_hazard)
dist_n_wantedS(1, last_hazard) = rest_sum
'Dim nmiss As Integer
'nmiss = 0
'Do Until n_actualS >= n_wantedS
'nmiss = nmiss + 1
'If nmiss > 100 Then
' Printdok "Not able to find enough immigrants in year " & year & "!"
' Exit Do
'End If
'Debug.Print "year= " & year & " nmiss= " & nmiss
'-----------------------------------------------------
'Step to allow households to be sampled randomly
ReDim slumpordningA(1 To m_hcount) As Double
Dim born_abroad_h As Integer
Dim j As Long
j = 0
For i = 1 To m_hcount
born_abroad_h = 1
If h_abroad(i) = 1 Then 'living abroad
inr = h_first_indnr(i)
Do Until inr = 0
index = indnr2index(inr)
' If all adult household members are born abroad the household
' is considered born abroad
If i_bvux(index) = 1 Then
born_abroad_h = mini(born_abroad_h, i_born_abroad(index))
End If
inr = i_next_indnr(index)
Loop
If born_abroad_h = 0 Then 'Born in Sweden
j = j + 1
slumpordningA(j) = i
End If
End If
Next
ReDim Preserve slumpordningA(1 To j)
ReDim slumpordning(1 To j, 1 To 2) As Single
For i = 1 To j
slumpordning(i, 1) = Rnd()
slumpordning(i, 2) = slumpordningA(i)
Next
'Call SORTMATRIX(j, 2, slumpordning(1, 1))
If j <= 100000 Then
Call SORTMATRIX(j, 2, slumpordning(1, 1))
Else
Call quicksort(slumpordning, 1)
End If
For hhi = 1 To j 'm_hcount
hh = slumpordning(hhi, 2)
tid = year - h_emig_year(hh)
If h_abroad(hh) = 1 And tid > 0 Then
rand = Rnd()
' tid = year - h_emig_year(hh)
If tid > last_hazard Then tid = last_hazard 'No hazard for times > 100 years
' If tid > 100 Then tid = 100 'No hazard for times > 100 years
'Emigrationsmodulen körs först, det innebär att individer inte kan emigrera samma år som
'de immigrerat. Eftersom hazard(0)=0 betyder det också att de inte kan immigrera samma år
'som de emigrerat.
If dist_n_wantedS(2, tid) <= dist_n_wantedS(1, tid) Then
'If rand < hazard(tid) Then
h_abroad(hh) = 0
h_new_immig(hh) = 1
n_actualS = n_actualS + h_size(hh) 'summera hushållets storlek
dist_n_wantedS(2, tid) = dist_n_wantedS(2, tid) + h_size(hh)
h_new_housing(hh) = 1 ' Marked for tenure choice
'loop över individerna i hushållet
inr = h_first_indnr(hh)
Do Until inr = 0
index = indnr2index(inr)
i_abroad(index) = 0
i_new_immig(index) = 1
i_binvar(index) = base_year + model_time
inr = i_next_indnr(index)
Loop
'-- OS050518: Kommenterar in gammalt villkor
If n_actualS >= n_wantedS Then Exit For
End If
End If
Next
'Debug.Print "year= " & year & " Immig Sv Diffen blev " & n_actualS - n_wantedS & " " & n_wantedS
'Loop
'Det följande steget markarer vilka hushåll som uppfyller villkoret för att kunna
'klonas genom att tilldela mark_h(i) värdet 1.
'If n_wantedI >= 1 Then
For i = 1 To m_hcount
mark_h(i) = 0
Next
Dim nmiss As Long
nmiss = 0
'---------------------------------------------------------------------
'Här specificeras fördelningen av invandrarfamiljer. Kolumnerna 1 och 2
'anger åldersintervall, kolumn 3 hushållsstorlek och kolumn fyra frekvensen
'av kombinationen ålder och hushållsstorlek. Ålder är hushållets max-ålder.
'Den absolut bästa lösningen för kloningen av invandrarfamiljer är att att
'skapa en pool av invandrare att slumpmässigt dra personer ur. Hur gör man det på bästa
'sätt? Går det exempelvis att göra detta som för emigranturvalet? Tanken är att lägga till
'Y personer till grunddata som markeras på ett sätt att de inte räknas under simuleringarna
'men att de alltid kan klonas liksom vilken annan familj som helst? De skulle även kunna
'ligga i två nya databaser, ii_inv.bin och hh_inv.bin, vilka vi kan kalla på vid
'kloningen istället. Strukturen på dessa skulle vara exakt desamma som för ii.bin och
'hh.bin.
'Dim invandring(1 To 17, 1 To 4) As Double
'
'invandring(1, 1) = 18: invandring(1, 2) = 24: invandring(1, 3) = 1: invandring(1, 4) = 0.158836689
'invandring(2, 1) = 25: invandring(2, 2) = 34: invandring(2, 3) = 1: invandring(2, 4) = 0.082774049
'invandring(3, 1) = 25: invandring(3, 2) = 34: invandring(3, 3) = 2: invandring(3, 4) = 0.055928412
'invandring(4, 1) = 25: invandring(4, 2) = 34: invandring(4, 3) = 3: invandring(4, 4) = 0.049217002
'invandring(5, 1) = 25: invandring(5, 2) = 34: invandring(5, 3) = 4: invandring(5, 4) = 0.026845638
'invandring(6, 1) = 35: invandring(6, 2) = 44: invandring(6, 3) = 1: invandring(6, 4) = 0.060402685
'invandring(7, 1) = 35: invandring(7, 2) = 44: invandring(7, 3) = 2: invandring(7, 4) = 0.033557047
'invandring(8, 1) = 35: invandring(8, 2) = 44: invandring(8, 3) = 3: invandring(8, 4) = 0.035794183
'invandring(9, 1) = 35: invandring(9, 2) = 44: invandring(9, 3) = 4: invandring(9, 4) = 0.049217002
'invandring(10, 1) = 45: invandring(10, 2) = 54: invandring(10, 3) = 1: invandring(10, 4) = 0.080536913
'invandring(11, 1) = 45: invandring(11, 2) = 54: invandring(11, 3) = 2: invandring(11, 4) = 0.044742729
'invandring(12, 1) = 45: invandring(12, 2) = 54: invandring(12, 3) = 3: invandring(12, 4) = 0.024608501
'invandring(13, 1) = 45: invandring(13, 2) = 54: invandring(13, 3) = 4: invandring(13, 4) = 0.031319911
'invandring(14, 1) = 55: invandring(14, 2) = 64: invandring(14, 3) = 1: invandring(14, 4) = 0.085011186
'invandring(15, 1) = 55: invandring(15, 2) = 64: invandring(15, 3) = 2: invandring(15, 4) = 0.044742729
'invandring(16, 1) = 65: invandring(16, 2) = 150: invandring(16, 3) = 1: invandring(16, 4) = 0.111856823
'invandring(17, 1) = 65: invandring(17, 2) = 150: invandring(17, 3) = 2: invandring(17, 4) = 0.024608501
' Adjustment of wanted number of first time immigrants due to discrepancies between
' wanted and realized number of swedish immigrants
' Note that the matrix of wanted firt time immigrants by age and family type
' also is proportionally adjusted
Dim reldiff_wantedI As Double
reldiff_wantedI = ((n_wantedS + n_wantedI) - n_actualS) / n_wantedI
n_wantedI = (n_wantedS + n_wantedI) - n_actualS
ReDim array_wantedI(1 To UBound(n_im_trans, 2), 1 To 2)
For i = 1 To UBound(n_im_trans, 2)
array_wantedI(i, 1) = n_wantedI * n_im_trans(4, i) * reldiff_wantedI
array_wantedI(i, 2) = 0
Next
ReDim slumpordningA(1 To m_hcount) As Double
'Dim j As Long
j = 0
For i = 1 To m_hcount
If h_abroad(i) = 0 And mark_h(i) = 0 Then 'And invandrings_year(h_hhnr(i)) > 0 Then
j = j + 1
slumpordningA(j) = i
End If
Next
ReDim Preserve slumpordningA(1 To j)
ReDim slumpordning(1 To j, 1 To 2) As Single
For i = 1 To j
slumpordning(i, 1) = Rnd()
slumpordning(i, 2) = slumpordningA(i)
Next
'Call SORTMATRIX(j, 2, slumpordning(1, 1))
If j <= 100000 Then
Call SORTMATRIX(j, 2, slumpordning(1, 1))
Else
Call quicksort(slumpordning, 1)
End If
Dim hushållsstorlek As Integer
Dim maxage As Integer
Dim htyp As Byte
n_actualI = 0
For hhi = 1 To j 'm_hcount
hh = slumpordning(hhi, 2)
'Följande är en (ful) kodslinga för att välja ut önskat antal familjer
'som uppfyller max_ålder och hushållsstorlek för kloning.
hushållsstorlek = h_size(hh)
'-- OS050517: Ändrade villkor för att bryta loopen
If Abs(n_wantedI - (n_actualI + hushållsstorlek)) > n_wantedI - n_actualI Then Exit For
maxage = max_age(h_hhnr(hh))
If hushållsstorlek = 1 And maxage > 17 And maxage < 25 Then
htyp = 1
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek = 2 And maxage > 17 And maxage < 25 Then
htyp = 2
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek = 3 And maxage > 17 And maxage < 25 Then
htyp = 3
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek >= 4 And maxage > 17 And maxage < 25 Then
htyp = 4
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek = 1 And maxage > 24 And maxage < 35 Then
htyp = 5
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek = 2 And maxage > 24 And maxage < 35 Then
htyp = 6
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek = 3 And maxage > 24 And maxage < 35 Then
htyp = 7
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek >= 4 And maxage > 24 And maxage < 35 Then
htyp = 8
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek = 1 And maxage > 34 And maxage < 45 Then
htyp = 9
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek = 2 And maxage > 34 And maxage < 45 Then
htyp = 10
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek = 3 And maxage > 34 And maxage < 45 Then
htyp = 11
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek >= 4 And maxage > 34 And maxage < 45 Then
htyp = 12
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek = 1 And maxage > 44 And maxage < 55 Then
htyp = 13
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek = 2 And maxage > 44 And maxage < 55 Then
htyp = 14
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek = 3 And maxage > 44 And maxage < 55 Then
htyp = 15
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek >= 4 And maxage > 44 And maxage < 55 Then
htyp = 16
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek = 1 And maxage > 54 And maxage < 65 Then
htyp = 17
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek = 2 And maxage > 54 And maxage < 65 Then
htyp = 18
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek = 3 And maxage > 54 And maxage < 65 Then
htyp = 19
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek >= 4 And maxage > 54 And maxage < 65 Then
htyp = 20
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek = 1 And maxage > 64 Then
htyp = 21
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek = 2 And maxage > 64 Then
htyp = 22
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek = 3 And maxage > 64 Then
htyp = 23
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
ElseIf hushållsstorlek >= 4 And maxage > 64 Then
htyp = 24
If array_wantedI(htyp, 1) > array_wantedI(htyp, 2) Then
array_wantedI(htyp, 2) = array_wantedI(htyp, 2) + hushållsstorlek
mark_h(hh) = 1
n_actualI = n_actualI + hushållsstorlek
End If
Else
End If
'*** Gammalt villkor If n_actualI >= n_wantedI Then Exit For
Next
' Clone the selected housholds
Call clone_hh
'Debug.Print "year= " & year & " Diffen blev " & n_actualI - n_wantedI
status "Immigration: " & (n_actualS + n_actualI) * m_weight & " (" _
& n_im(mini(2110, year), 1) + n_im(mini(2110, year), 2) & " )"
m_immigrated = (n_actualS + n_actualI) * m_weight
End Sub
Public Sub emigration()
'I denna rutin ska jag göra följande:
' 1) Antag bruttoemigration, Y, fördelat på Yi och Ys, vilka
' hämtas ifrån ett Excel-ark
' 2) Simulera vilka familjer som lämnar Sverige. Avgör också vilka familjer
' som vi måste hålla reda på för dels framtida återinvandring, dels
' för pensionsutbetalningar utomlands.
Dim n_wantedS As Long
Dim n_wantedI As Long
Dim n_actualS As Long
Dim n_actualI As Long
Dim year As Integer
Dim hh As Long
Dim hhi As Long
'Dim inr As Long
'Dim index As Long
Dim rand As Double
Dim pred As Double
Dim scale_par As Double
'Alignment data constant from 2050
year = base_year + model_time
scale_par = 1 / m_weight
'n_wantedS = scale_par * n_em(year + 1, 1) 'Migrationsdata börjar på 1998, är det ett problem???
n_wantedS = scale_par * n_em(mini(2110, year), 1) 'Migrationsdata börjar på 1998, är det ett problem???
n_wantedI = scale_par * n_em(mini(2110, year), 2) 'Migrationsdata börjar på 1998, är det ett problem???
n_actualS = 0 'n_actual personer som emigrerat i snurran nedan.
n_actualI = 0 'n_actual personer som emigrerat i snurran nedan.
Dim tid_i_sverige As Integer
Dim age_max As Integer
Dim edu_max As Integer
Dim binvar_max As Integer
Dim vuxna As Integer
Dim barn As Integer
Dim e2 As Integer
Dim e3 As Integer
Dim ex As Double
Dim nmiss As Integer
nmiss = 0
Do Until n_actualS >= n_wantedS And n_actualI >= n_wantedI
nmiss = nmiss + 1
If nmiss > 100 Then
If nmiss = 101 Then nmiss = 100
Printdok "Not able to find enough emigrants! Please consult programmer."
Exit Do
End If
'-----------------------------------------------------
'Step to allow households to be sampled randomly
ReDim slumpordningA(1 To m_hcount) As Double
Dim j As Long
j = 0
For i = 1 To m_hcount
If h_abroad(i) = 0 Then
j = j + 1
slumpordningA(j) = i
End If
Next
ReDim Preserve slumpordningA(1 To j)
ReDim slumpordning(1 To j, 1 To 2) As Single
For i = 1 To j
slumpordning(i, 1) = Rnd()
slumpordning(i, 2) = slumpordningA(i)
Next
'Call SORTMATRIX(j, 2, slumpordning(1, 1))
If j <= 100000 Then
Call SORTMATRIX(j, 2, slumpordning(1, 1))
Else
Call quicksort(slumpordning, 1)
End If
'-----------------------------------------------------
Dim go As Integer 'Har lagts till för att garanatera att rätt person jämförs vid rand 0 And _
Abs(n_wantedI - (n_actualI + h_size(hh))) < n_wantedI - n_actualI Then
ex = Exp(-1.5157 + age_max * (-0.00774) + barn * (-0.5664) + _
vuxna * (-0.4291) + barn * vuxna * 0.2024 + e2 * (-0.7397) + _
e3 * (-0.542) + tid_i_sverige * (-0.0685))
pred = ex / (1 + ex)
go = 1
If rand < pred And go = 1 Then
Call emigrate_hh(hh, year)
n_actualI = n_actualI + h_size(hh)
End If
End If
If tid_i_sverige = -999 And _
Abs(n_wantedS - (n_actualS + h_size(hh))) < n_wantedS - n_actualS Then
ex = Exp(-4.5667 + age_max * (-0.0212) + barn * 0.0365 + _
vuxna * (-0.927) + barn * vuxna * 0.00988 + e2 * 0.278 + _
e3 * 1.2092)
pred = ex / (1 + ex)
go = 1
If rand < pred And go = 1 Then
Call emigrate_hh(hh, year)
n_actualS = n_actualS + h_size(hh)
End If
End If
' If rand < pred And go = 1 Then
'
' h_abroad(hh) = 1
' h_emig_year(hh) = year
' h_new_em(hh) = 1
' If tid_i_sverige > 0 And n_actualI < n_wantedI Then
' n_actualI = n_actualI + h_size(hh)
' End If
' If tid_i_sverige = -999 And n_actualS < n_wantedS Then
' n_actualS = n_actualS + h_size(hh)
' End If
'
' 'loop over individuals in household
' inr = h_first_indnr(hh)
' Do Until inr = 0
' index = indnr2index(inr)
' i_abroad(index) = 1
' i_emig_year(index) = year
' i_new_em(index) = 1
' inr = i_next_indnr(index)
' Loop
' End If
End If
Next
Loop
status "Emigration: " & (n_actualS + n_actualI) * m_weight & " (" & n_em(mini(2110, year), 1) + n_em(mini(2110, year), 2) & " )"
m_emigrated = (n_actualS + n_actualI) * m_weight
'Debug.Print "year= " & year & " S:" & n_actualS - n_wantedS & " I:" & n_actualI - n_wantedI
' Debug.Print "year= " & year & " Emigration= " & n_actualS '-m_netmigration
End Sub
Public Sub emigrate_hh(hh As Long, year As Integer)
Dim inr As Long, index As Long
h_abroad(hh) = 1
h_emig_year(hh) = year
' h_new_em(hh) = 1
'-- loop over individuals in household
inr = h_first_indnr(hh)
Do Until inr = 0
index = indnr2index(inr)
i_abroad(index) = 1
i_emig_year(index) = year
i_new_em(index) = 1
inr = i_next_indnr(index)
Loop
End Sub
'Procedur för att plocka ut familjens maxvärden till skattningarna
'av emigrationssannolikhet. hhnr motsvarar hushållets nummer
Public Sub hh_mix(hhnr As Long, ByRef age_max As Integer, ByRef edu_max As Integer, _
ByRef vuxna As Integer, ByRef barn As Integer, _
ByRef binvar_max As Integer, ByRef tid_i_sverige As Integer)
Dim inr As Long
Dim index
Dim age As Integer
Dim edu As Integer
Dim binvar As Integer
Dim born_abroad_h As Integer
Dim born_abroad_h_max As Integer
vuxna = h_n_adults(hhnr2index(hhnr)) 'antal_vuxna
barn = h_n_child(hhnr2index(hhnr)) 'antal_barn
inr = h_first_indnr(hhnr2index(hhnr))
age_max = 0
edu_max = 0
binvar_max = 0
born_abroad_h_max = 0
Do Until inr = 0
index = indnr2index(inr)
age = i_age(index)
edu = i_edlevel(index)
binvar = i_binvar(index)
born_abroad_h = i_born_abroad(index)
If age > age_max Then age_max = age
If edu > edu_max Then edu_max = edu
If binvar > binvar_max Then binvar_max = binvar
If born_abroad_h > born_abroad_h_max Then born_abroad_h_max = born_abroad_h
inr = i_next_indnr(index)
Loop
If born_abroad_h_max = 1 And binvar_max > 0 Then 'Testa om det finns personer med binvar>0 och born_abroad=1
tid_i_sverige = base_year + model_time - binvar_max
Else
tid_i_sverige = -999
End If
End Sub
'Funktion för att plocka invandringsår för hushåll I. Vi antar att högsta året
'gäller för hela familjen. hhnr är hushållsnumret
Public Function invandrings_year(hhnr As Long) As Integer
Dim inr As Long
Dim index As Long
Dim binvar As Integer
Dim binvar_max As Integer
inr = h_first_indnr(hhnr2index(hhnr))
binvar_max = 0
Do Until inr = 0
index = indnr2index(inr)
binvar = i_binvar(index)
If binvar > binvar_max Then binvar_max = binvar
inr = i_next_indnr(index)
Loop
invandrings_year = binvar_max
End Function
'Function to evaluate age of oldest individual in family hhnr
Public Function max_age(hhnr As Long) As Integer
Dim inr As Long
Dim index As Long
inr = h_first_indnr(hhnr2index(hhnr))
max_age = 0
Do Until inr = 0
index = indnr2index(inr)
max_age = maxi(max_age, i_age(index))
inr = i_next_indnr(index)
Loop
End Function
'**********************************************************************************
' Sub cohab selects single adult females (not living abroad) to form new households.
' The selected women are matched against all available men (single, adult and
' not living abroad) to find the most suitable man (smallest age-difference).
' Variance reduction is applied so that the number of new cohabitations (= # of
' selected women) equals the expected value of the process.
'**********************************************************************************
Public Sub cohab()
'! Handles cohabitation
Printdok "cohab"
Dim i As Long, at_risk As Long
Dim prob As Double, expval As Double, rand() As Double
Dim sortvek() As Double, probs() As Double, xbeta As Double
Dim pool_female_atrisk() As Long
Dim age As Integer
'*** Special calculation for validation of pension system rules
'*** TP030212
If get_scalefactor_active("Pension_Orange") = 1 Then
status "Orange: no cohab!"
Else
n = 0
' Create pool of available men
' NOTE: Adult, living in single households in Sweden (no emigrants)
pool_male_size = 0
For i = 1 To m_icount
If i_sex(i) = 1 And i_bvux(i) = 1 And _
h_n_adults(hhnr2index(i_hhnr(i))) = 1 And i_abroad(i) = 0 Then
pool_male_size = pool_male_size + 1
pool_male(pool_male_size) = i_indnr(i)
End If
Next
' 1) create a pool of women at risk and calculate the expected number of selected women
expval = 0
at_risk = 0
ReDim pool_female_atrisk(1 To m_icount), probs(1 To m_icount)
For i = 1 To m_icount
' Population at risk: women living alone in own household in Sweden (no emigrants)
If i_sex(i) = 2 And i_bvux(i) = 1 And _
h_n_adults(hhnr2index(i_hhnr(i))) = 1 And i_abroad(i) = 0 Then
at_risk = at_risk + 1
pool_female_atrisk(at_risk) = i_indnr(i)
'*** The probability that the woman enters a cohabitation is modeled using
'*** a logistic regression with age as only explanatory variable.
'*** NOTE: for purposes of alignment to observed data individuals aged 23 and
'*** younger have probabilities equal to individuals aged 24. Individuals aged
'*** 61 and above have the same probabilities as individuals aged 60. (TP 041028)
age = mini(60, maxi(21, i_age(i)))
xbeta = -41.94824142 + 0.4 + _
(age / 10) * 49.39597128 + _
(age / 10) ^ 2 * -21.85389201 + _
(age / 10) ^ 3 * 4.086284708 + _
(age / 10) ^ 4 * -0.27641113
prob = 1 / (1 + Exp(-xbeta))
'*** The probability cannot be exactly zero or one because of the logit
'*** transformation below
probs(at_risk) = maxi(0.00001, mini(0.99999, prob))
expval = expval + probs(at_risk) * m_weight
End If
Next
' Delete empty space
ReDim Preserve pool_female_atrisk(1 To at_risk)
ReDim Preserve probs(1 To at_risk)
' 2) Draw uniform numbers
ReDim rand(1 To at_risk)
Call RANUNI(at_risk, rand(1), model_time + base_year + random * Rnd)
' 3) Calculate the logit transformed sorting variable
' NOTE: negative sign in order to get the right selection in findlimit
ReDim sortvek(1 To at_risk)
For i = 1 To at_risk
sortvek(i) = -(logit(rand(i)) - logit(probs(i)))
Next
' 4) Find the threshold value for choosing the right number of women
' NOTE: naive simulation uses limit = 0!!
Dim limit As Double
limit = findlimit(sortvek, CLng(expval / m_weight), 0)
' 5) Create pool of selected women
pool_female_size = 0
' Dim fem_ages() As Integer
' ReDim fem_ages(1 To m_icount)
For i = 1 To at_risk
If sortvek(i) > limit Then
pool_female_size = pool_female_size + 1
pool_female(pool_female_size) = pool_female_atrisk(i)
' fem_ages(pool_female_size) = i_age(indnr2index(pool_female_atrisk(i)))
End If
Next
' MsgBox PrintFile("c:\slask.txt", fem_ages)
status "Cohab " & pool_female_size * m_weight
If pool_female_size > 0 Then Call match_couples
End If
End Sub
'******************************************************************************************
'*** Sub early_retire randomizes individuals into disability. Complete documentation of the
'*** estimated models can be found at G:\Sesim\Dokumentation\Estimering\Disability pension.
'*** The predicted number of newly disabled is aligned to actual data (se reference below)
'*** until 2001 and thereafter to the expected number of newly disabled.
'******************************************************************************************
Public Sub early_retire()
'! Routine for early retirement
Printdok "early_retire"
Dim i As Long, n As Long
Dim probs() As Double, rand() As Double, sortvek() As Double
Dim pgi_atrisk() As Double, pool_atrisk() As Long
ReDim probs(1 To m_icount), pool_atrisk(1 To m_icount), rand(1 To m_icount)
ReDim sortvek(1 To m_icount), pgi_atrisk(1 To m_icount)
Dim count_atrisk As Long, count_er As Long
Dim expected As Double
count_atrisk = 0
expected = 0
count_er = 0
'*** Step 1: create pool of individuals at risk and calculate quartiles of PGI
For i = 1 To m_icount
'*** Individuals at risk
If i_age(i) >= 16 And i_age(i) <= 64 And i_status(i) <> 4 And i_abroad(i) = 0 Then
count_atrisk = count_atrisk + 1
pool_atrisk(count_atrisk) = i_indnr(i) '*** add individuals at risk to pool
pgi_atrisk(count_atrisk) = CDbl(i_pgi(i)) '*** add PGI to pool
End If
Next
'*** Pack and/or define vectors depending on calculated count_atrisk
ReDim Preserve pool_atrisk(1 To count_atrisk), pgi_atrisk(1 To count_atrisk)
ReDim Preserve probs(1 To count_atrisk), rand(1 To count_atrisk)
ReDim Preserve sortvek(1 To count_atrisk)
'*** Calculate limits for quartiles of PGI for individuals at risk
Dim pgi_q() As Double
pgi_q = arr_Percentile(pgi_atrisk, 25, 50, 75)
Dim idx As Long
Dim xbeta As Double
'*** Step 2: calculate individual probabilities for entering disability
Dim ed16_29, ed30_60, pgi16_29, pgi30_60, pgi_bkon_30_60, bald_pgi
Dim pgi61_, ed61_, pgi_bkon_61_, ed_bkon_61_, bald_ed_61_
Dim bald16 As Double, icp16_29 As Double, icp30_60 As Double, icp_61_ As Double
Dim bald16_29 As Double, bald30_60 As Double, bald61_ As Double, bkon61 As Double
Dim bald61 As Double, icp61_ As Double
Dim bald16_ As Integer, pgi_ As Integer, bkon_ As Integer, civstat_ As Integer, inv_ As Integer
Dim bkon30_60 As Double, civstat As Double, inv As Double, bald_civstat As Double
'*** Parameter estimates
'*** 16-29
icp16_29 = -12.20176621
ed16_29 = Array(0.953595481, -0.069001588, -0.884593893)
bald16_29 = 0.177644219
bald16 = 1.89992357
pgi16_29 = Array(1.627132586, 0.547965417, -0.286313535, -1.888784468)
'*** 30-60
icp30_60 = -10.30358891
ed30_60 = Array(0.214902421, 0.105155022, -0.320057443)
pgi30_60 = Array(3.237931598, 1.928913217, -1.036837021, -4.130007795)
bkon30_60 = -0.006222867
pgi_bkon_30_60 = Array(-0.286842131, 0.097124963, 0.17195553, 0.017761638) '*** ref: men
civstat = 0.908264782
inv = -0.150018962
bald30_60 = 0.110462899
bald_pgi = Array(-0.04938895, -0.022603466, 0.015495372, 0.056497044)
bald_civstat = -0.014096414 '*** ref: non-cohab
'*** 60-
icp61_ = 10.4165131
pgi61_ = Array(-0.431271281, 0.558623526, 0.301083363, -0.428435607)
ed61_ = Array(-12.36657816, -4.49624942, 16.86282758)
bkon61 = -0.037840852
pgi_bkon_61_ = Array(-0.405932962, -0.02343662, 0.178264996, 0.251104586)
ed_bkon_61_ = Array(0.244462498, -0.026095518, -0.218366981)
bald61 = -0.229034907
bald_ed_61_ = Array(0.201916564, 0.071616875, -0.27353344)
For i = 1 To count_atrisk
idx = indnr2index(pool_atrisk(i))
If i_pgi(idx) < pgi_q(1, 2) Then
pgi_ = 1
Else
If i_pgi(idx) < pgi_q(2, 2) Then
pgi_ = 2
Else
If i_pgi(idx) < pgi_q(3, 2) Then
pgi_ = 3
Else
pgi_ = 4
End If
End If
End If
If i_sex(idx) = 1 Then bkon_ = 1 Else bkon_ = -1
If i_civ_stat(idx) = 0 Then civstat_ = 1 Else civstat_ = -1
If i_born_abroad(idx) = 0 Then inv_ = 1 Else inv_ = -1
Select Case i_age(idx)
'*** 16 - 29 years of age
Case Is <= 29
If i_age(idx) = 16 Then bald16_ = 1 Else bald16_ = 0
xbeta = icp16_29 + ed16_29(i_edlevel(idx)) + i_age(idx) * bald16_29 + _
bald16 * bald16_ + pgi16_29(pgi_ - 1)
'*** 30 - 60 years of age
Case 30 To 60
xbeta = icp30_60 + pgi30_60(pgi_ - 1) + ed30_60(i_edlevel(idx)) + _
bkon_ * bkon30_60 + pgi_bkon_30_60(pgi_ - 1) * bkon_ + civstat_ * civstat + _
inv_ * inv + i_age(idx) * bald30_60 + bald_pgi(pgi_ - 1) * i_age(idx) + _
bald_civstat * i_age(idx) * civstat_
'*** 61 - years of age
Case Is >= 61
xbeta = icp61_ + pgi61_(pgi_ - 1) + ed61_(i_edlevel(idx)) + bkon61 * bkon_ + _
pgi_bkon_61_(pgi_ - 1) * bkon_ + ed_bkon_61_(i_edlevel(idx)) * bkon_ + _
bald61 * i_age(idx) + bald_ed_61_(i_edlevel(idx)) * i_age(idx)
End Select
' Calculate probability
probs(i) = 1 / (1 + Exp(-xbeta))
' Ackumulate expected count of newly disabled
expected = expected + probs(i)
Next
'*** Step 3: draw random numbers
Call RANUNI(count_atrisk, rand(1), model_time + base_year + random * Rnd)
'*** Step 4: calculate the logit transformed "sorting" variable
'*** NOTE: negative sign in order to get the right selection in findlimit
For i = 1 To count_atrisk
sortvek(i) = -(logit(rand(i)) - logit(probs(i)))
Next
'*** From base_year until the current year there are real data to align the number of
'*** new disability pensioneers to. (Data source: RFV, statistikinformation Is-I 2002:1,
'*** tabell 1 - Nybeviljade förtidspensioner och sjukbidrag med fördelning efter
'*** pensioneringsgrad och kön.)
'*** 2002-2008: Nytillkomna (i sjukers-de från aktivers) enligt budgetunderlag korrigerat för heltid enl
' (Helårsekvivalenter 2004/Antal 2004 enl budgetunderlag) = 0.84
'*** 060207 ThP Kompleterad med utfall t.o.m. 2005 se;
'*** S:\Dokument\Dokumentation\SESIM\Alignment\Nybeviljade SA per ålder från FK.xls
Select Case (base_year + model_time)
Case 2000
expected = 40845 / m_weight
Case 2001
expected = 47501 / m_weight
Case 2002
expected = 52648 / m_weight
Case 2003
expected = 52499 / m_weight
Case 2004
expected = 56916 / m_weight
Case 2005
expected = 46548 / m_weight
Case 2006
expected = 37550 / m_weight
Case Else
End Select
'*** Step 5: Find the threshold value for choosing the right number of women
'*** NOTE: naive simulation uses limit = 0!
Dim limit As Double
limit = findlimit(sortvek, CLng(expected), 0)
'*** Step 6: Loop through the pool of individuals at risk and select the ones that
'*** will be randomized into disability
For i = 1 To count_atrisk
i_new_fp(indnr2index(pool_atrisk(i))) = 0
If sortvek(i) > limit Then
i_new_fp(indnr2index(pool_atrisk(i))) = 1 '*** disabled
count_er = count_er + 1 '*** counts # of newly disabled
i_ftp_year(indnr2index(pool_atrisk(i))) = base_year + model_time
End If
Next
status "Early retired: " & count_er * m_weight
End Sub
'****************************************************************************************
'*** Sub rehab randomizes individuals to exit the disability status (rehabilitation).
'*** Complete documentation of the estimated models can be found at
'*** G:\Sesim\Dokumentation\Estimering\Disability pension.
'*** Since no reliable data on the number of rehabilitations have been found the model
'*** aligns the number of rehabilitated individuals to the expected number.
'****************************************************************************************
Public Sub rehab()
'! Routine for rehabilition from early retirement
Printdok "rehab"
Dim expected As Double, xbeta As Double, limit As Double
Dim count_atrisk As Long, count_event As Long, i As Long, duration As Long, idx As Long
Dim probs() As Double, rand() As Double, sortvek() As Double
Dim pool_atrisk() As Long
ReDim probs(1 To m_icount), pool_atrisk(1 To m_icount)
expected = 0
count_atrisk = 0
count_event = 0
Dim ed
ed = Array(-0.616181648, 0.101532353, 0.514649295) ' edlevel parameter
'*** Step 1: create pool of individuals at risk, calculate probabilities and expected number
'*** of events.
'*** NOTE: only individuals that were disabled at the end of last year is at risk. Hence, no
'*** one can be disabled and rehabilitated during one year.
For i = 1 To m_icount
If 16 <= i_age(i) And i_age(i) <= 64 And i_status(i) = 4 Then
count_atrisk = count_atrisk + 1
pool_atrisk(count_atrisk) = i_indnr(i)
' duration truncated at 10 years
'*** OBSERVERA: i_ftp_year anges med två siffror nu - blir fyra siffror i nästa version
'*** av startdata!!!
If i_ftp_year(i) < 100 Then
duration = mini(10, ((base_year + model_time) - (1900 + i_ftp_year(i))))
Else
duration = mini(10, ((base_year + model_time) - i_ftp_year(i)))
End If
xbeta = -0.988939296 - 0.148911411 * duration - 0.079384705 * i_age(i) + ed(i_edlevel(i))
probs(count_atrisk) = 1 / (1 + Exp(-xbeta))
expected = expected + probs(count_atrisk)
End If
Next
'*** Pack vectors
ReDim Preserve pool_atrisk(1 To count_atrisk), probs(1 To count_atrisk)
ReDim sortvek(1 To count_atrisk), rand(1 To count_atrisk)
'*** Step 2: draw uniform random numbers
Call RANUNI(count_atrisk, rand(1), model_time + base_year + random * Rnd)
'*** Step 3: calculate the logit transformed "sorting" variable
'*** NOTE: negative sign in order to get the right selection in function findlimit
For i = 1 To count_atrisk
sortvek(i) = -(logit(rand(i)) - logit(probs(i)))
Next
'*** Step 4: find the threshold value for choosing the right number of events
limit = findlimit(sortvek, CLng(expected), 0)
'*** Step 5: loop through the pool of individuals at risk and select the ones that
'*** will be rehabilitated from disability
For i = 1 To count_atrisk
If sortvek(i) > limit Then
count_event = count_event + 1
idx = indnr2index(pool_atrisk(i))
i_new_fp(idx) = -1
End If
Next
status "Rehab: " & count_event * m_weight
End Sub
Public Sub demograf_stat()
'! Updates some variables after demographics
Printdok "demograf_stat"
Dim indnr As Long, idnr_h As Long, index_h As Long, idnr_i As Long, index_i As Long
Dim pop1(0 To 106, 1 To 2) As Double, pop2(0 To 106, 1 To 2) As Double
'*** Random numbers for assignment of permanent "luck factors" for
'*** statistical panel data models
Dim randnr() As Double
Dim randcount As Long, drawcount As Integer
ReDim randnr(1 To m_icount)
Call RANNOR(m_icount, randnr(1), model_time + base_year + random * Rnd)
randcount = 1
drawcount = 1
'*** Loop across households and individuals
For index_h = 1 To m_hcount
h_max_age(index_h) = 0
h_max_edlevel(index_h) = 0
h_bvux_work(index_h) = 1
h_max_inc_taxable(index_h) = 0
h_indnr_female(index_h) = 0
h_indnr_male(index_h) = 0
h_n_childlt7(index_h) = 0
h_n_childlt12(index_h) = 0
idnr_i = h_first_indnr(index_h)
index_i = indnr2index(idnr_i)
h_n_ge18(index_h) = 0
Do While idnr_i <> 0
index_i = indnr2index(idnr_i)
'*** Calculate maximum age in household
If i_age(index_i) > h_max_age(index_h) Then
h_max_age(index_h) = i_age(index_i)
End If
'*** Calculate maximum educational attainment in household
If i_edlevel(index_i) > h_max_edlevel(index_h) Then
h_max_edlevel(index_h) = i_edlevel(index_i)
End If
If i_bvux(index_i) = 1 And i_status(index_i) <> 8 Then
h_bvux_work(index_h) = 0
End If
'*** Calculate maximum taxable income in household
If i_inc_taxable1(index_i) > h_max_inc_taxable(index_h) Then
h_max_inc_taxable(index_h) = i_inc_taxable1(index_i)
End If
'*** Calculate individual numbers of adult male and female in household
If i_bvux(index_i) = 1 And i_sex(index_i) = 1 Then
h_indnr_male(index_h) = i_indnr(index_i)
End If
If i_bvux(index_i) = 1 And i_sex(index_i) = 2 Then
h_indnr_female(index_h) = i_indnr(index_i)
End If
'*** Accumulate number of children under the age of seven
If i_age(index_i) < 7 Then
h_n_childlt7(index_h) = h_n_childlt7(index_h) + 1
End If
If i_age(index_i) < 12 Then
h_n_childlt12(index_h) = h_n_childlt12(index_h) + 1
End If
If i_abroad(index_i) = 0 Then i_botid(index_i) = i_botid(index_i) + 1
'*** counts of individuals in Sweden and abroad by sex and age (0 - 106+)
If i_abroad(index_i) = 0 Then
pop1(mini(106, i_age(index_i)), i_sex(index_i)) = _
pop1(mini(106, i_age(index_i)), i_sex(index_i)) + m_weight
Else
pop2(mini(106, i_age(index_i)), i_sex(index_i)) = _
pop2(mini(106, i_age(index_i)), i_sex(index_i)) + m_weight
End If
If i_age(index_i) >= 18 Then h_n_ge18(index_h) = h_n_ge18(index_h) + 1
'*** If all random numbers are used up - make new draw
If randcount = m_icount Then
Call RANNOR(m_icount, randnr(1), model_time + drawcount + random * Rnd)
drawcount = drawcount + 1
randcount = 1
End If
'*** Set individual random factors for panel data models
If i_unemp_ivariance(index_i) = 0 Then
i_unemp_ivariance(index_i) = randnr(randcount) * Sqr(0.648)
randcount = randcount + 1
End If
If i_inc_ivariance(index_i) = 0 Then
i_inc_ivariance(index_i) = randnr(randcount) * Sqr(0.0773938)
randcount = randcount + 1
End If
If i_sickleave_ivariance(index_i) = 0 Then
i_sickleave_ivariance(index_i) = randnr(randcount) * Sqr(0.638)
randcount = randcount + 1
End If
idnr_i = i_next_indnr(index_i)
Loop
' Determine index for head of household
Select Case h_n_adults(index_h)
Case 1
If h_indnr_female(index_h) <> 0 Then
h_idx_headofhh(index_h) = indnr2index(h_indnr_female(index_h))
Else
h_idx_headofhh(index_h) = indnr2index(h_indnr_male(index_h))
End If
Case 2
If i_age(indnr2index(h_indnr_female(index_h))) > _
i_age(indnr2index(h_indnr_male(index_h))) Then
h_idx_headofhh(index_h) = indnr2index(h_indnr_female(index_h))
Else
h_idx_headofhh(index_h) = indnr2index(h_indnr_male(index_h))
End If
Case Else
status "ERROR (demograf_stat): invalid number of adults!"
End Select
'*** Household size
h_size(index_h) = get_hh_size(h_hhnr(index_h))
'*** Nr of children
h_n_child(index_h) = h_size(index_h) - h_n_adults(index_h)
'*** Municipality code
h_kommunkod(index_h) = kommundata(h_kommunindex(index_h)).kod
' h_la_region(index_h) = kommundata(h_kommunindex(index_h)).LA_region
h_BB_region(index_h) = BB_Region(h_kommunkod(index_h))
Next
'*** write demographic information to textfile as default. This is easily imported to SAS
'*** or Excel for further analysis
Dim demofile As Integer
demofile = FreeFile
'*** Create the file at base-year and append information for each year in simulation
If model_time = 0 Then
Open sesimpath & "\tempdata\demo_count.txt" For Output As #demofile
Else
Open sesimpath & "\tempdata\demo_count.txt" For Append As #demofile
End If
If model_time = 0 Then _
Print #demofile, "year", "age", "men_swe", "women_swe", "men_abroad", "women_abroad"
For i = 0 To 106
Print #demofile, base_year + model_time, i, pop1(i, 1), pop1(i, 2), pop2(i, 1), pop2(i, 2)
Next i
Close #demofile
End Sub
'********************************************************************************
' Subroutine Split_Pair determines which couples that are to separate during the
' year.
' The Monte-Carlo variance is eliminated by aligning the total number of
' separations to the expected number of separations.
'********************************************************************************
Public Sub Split_Pair()
Printdok "split_pair"
Dim i As Long, n As Long, prob As Double, limit As Double
Dim expval As Double, probs() As Double, randvek() As Double, r() As Double
Dim xbeta As Double
Dim pool_atrisk() As Long, n_atrisk As Long, hhidx As Long
Dim age As Integer
'*** Special calculation for validation of pension system rules
'*** TP030212
If get_scalefactor_active("Pension_Orange") = 1 Then
status "Orange: no split_pair!"
Else
Call clear_marks '*** This needs to be checked again! TP 020207
expval = 0
n_atrisk = 0
ReDim pool_atrisk(1 To m_icount)
ReDim probs(1 To m_icount)
'*** 1) Create pool of individuals at risk, calculate the expected number of
'*** separations and the probability vector
For i = 1 To m_icount
hhidx = hhnr2index(i_hhnr(i))
'*** Population at risk: cohabiting men who didn't enter the cohabitation
'*** during the current year
If i_sex(i) = 1 And i_bvux(i) = 1 And h_n_adults(hhidx) = 2 And _
h_form_year(hhidx) < base_year + model_time Then
n_atrisk = n_atrisk + 1
pool_atrisk(n_atrisk) = i_indnr(i)
age = mini(i_age(i), 40) 'Constant risk from 40
xbeta = -326.6828906 + (age / 10) * 451.4279587 + _
(age / 10) ^ 2 * -230.9062246 + _
(age / 10) ^ 3 * 51.44270998 + _
(age / 10) ^ 4 * -4.222626013
prob = 1 / (1 + Exp(-xbeta))
' avoid problems with logit transformation later
probs(n_atrisk) = maxi(0.0000001, mini(0.9999999, prob))
' Ackumulate the expected number of separations
expval = expval + probs(n_atrisk)
End If
Next
' Pack vectors
ReDim Preserve probs(1 To n_atrisk)
ReDim Preserve pool_atrisk(1 To n_atrisk)
'*** 2) Draw random numbers
ReDim randvek(1 To n_atrisk)
ReDim r(1 To n_atrisk)
Call RANUNI(n_atrisk, randvek(1), model_time + base_year + random * Rnd)
'*** 3) Calculate the logit transformed sorting variable
' NOTE: reversed sign to find smallest values in function findlimit
For i = 1 To n_atrisk
r(i) = -(logit(randvek(i)) - logit(probs(i)))
Next
'*** 4) Find limit value for randomization
limit = findlimit(r, CLng(expval), 0)
'*** 5) Mark the randomized individuals/households
n = 0
Dim indidx As Long
For i = 1 To n_atrisk
If r(i) > limit Then
n = n + 1
indidx = indnr2index(pool_atrisk(i))
hhidx = hhnr2index(i_hhnr(indidx))
mark_i(indidx) = 1
mark_h(hhidx) = 1
End If
Next
status "Split pair " & n * m_weight
If n > 0 Then Call split_pair_doit
End If
End Sub
'***********************************************************************
' Sub split_pair_doit is breaking up the cohabitation for households
' selected in sub split_pair. Children stays with the women in the old
' household. Men are assigned to a new household
'***********************************************************************
Public Sub split_pair_doit()
Printdok "split_pair_doit"
Dim i As Long, j As Long, inr As Long, oldhnr As Long, inr_f As Long
Dim hhsize As Integer
Dim temp1 As Double, temp2 As Double
Dim oldhidx As Long
' Update pointers to next ind and pointers to first ind
For i = 1 To m_icount
' Men who are about to leave their households are marked
If mark_i(i) = 1 Then
oldhnr = i_hhnr(i)
oldhidx = hhnr2index(oldhnr)
' Check womens ind-nr (for lifeevent only)
inr = h_first_indnr(oldhidx)
inr_f = inr 'detta är en ny rad
Do Until i_next_indnr(indnr2index(inr)) = 0
If i_bvux(indnr2index(inr)) = 1 And i_sex(indnr2index(inr)) = 2 Then Exit Do
inr = i_next_indnr(indnr2index(inr))
inr_f = inr 'detta är också en ny rad
Loop
lifehist.write_hist i_indnr(i), "Man leaving women " & inr
lifehist.write_hist inr, "Man leaving: " & i_indnr(i)
' Updating civil status
i_civ_stat(indnr2index(inr)) = 0 ' female
i_civ_stat(i) = 0 ' male
' Year for separation
i_year_sep(indnr2index(inr)) = base_year + model_time
i_year_sep(i) = base_year + model_time
If (i_year_sep(indnr2index(inr)) = 1998 Or i_year_sep(i) = 1998) And _
base_year + model_time = 1999 Then
MsgBox "error... " & i
End If
' error handling
If h_form_year(hhnr2index(i_hhnr(i))) = 0 Then h_form_year(hhnr2index(i_hhnr(i))) = base_year
' Duration of household
i_last_hh_dur(indnr2index(inr)) = base_year + model_time - _
h_form_year(hhnr2index(i_hhnr(indnr2index(inr))))
i_last_hh_dur(i) = base_year + model_time - _
h_form_year(hhnr2index(i_hhnr(i)))
' Number of previous households / unions
i_nr_prev_hh(indnr2index(inr)) = i_nr_prev_hh(indnr2index(inr)) + 1
i_nr_prev_hh(i) = i_nr_prev_hh(i) + 1
hhsize = get_hh_size(oldhnr)
If hhsize = 1 Then
MsgBox "Error: Split hhsize=1! indnr=" & CStr(i_indnr(i)) & " hhnr=" & CStr(oldhnr)
Exit Sub
End If
' Disconnect the man from other individuals in household
'*** Case 1 - the man is first in the household
If i_indnr(i) = h_first_indnr(hhnr2index(oldhnr)) Then
' Change first_indnr in hh
h_first_indnr(hhnr2index(oldhnr)) = i_next_indnr(i)
'*** Case 2 - the man is last in the household
ElseIf i_next_indnr(i) = 0 Then
' change the zero pointer to the individual before the man
inr = h_first_indnr(hhnr2index(oldhnr))
For j = 2 To hhsize - 1
inr = i_next_indnr(indnr2index(inr))
Next
i_next_indnr(indnr2index(inr)) = 0
'*** Case 3 - the man is neither last nor first in the household
Else
' The individual before the man should point to the individual after the man
inr = h_first_indnr(hhnr2index(oldhnr))
Do Until i_next_indnr(indnr2index(inr)) = i_indnr(i)
inr = i_next_indnr(indnr2index(inr))
Loop
i_next_indnr(indnr2index(inr)) = i_next_indnr(i)
End If
'*** Update individual data for male
största_hhnr = största_hhnr + 1
i_hhnr(i) = största_hhnr
i_next_indnr(i) = 0
' Allocate more memory for household data if necessary
' NOTE: adding 500 elements at a time increases speed of computations
If m_hcount = UBound(h_hhnr) Then
Call dyn_vect_h(m_hcount + 500)
'Debug.Print "allocating memory - split_pair_doit"
End If
m_hcount = m_hcount + 1
'*** Copy household information from womens household and update some variables for the
'*** new single male household
h_hhnr(m_hcount) = största_hhnr
h_first_indnr(m_hcount) = i_indnr(i)
h_kommunindex(m_hcount) = h_kommunindex(oldhidx)
h_kommunkod(m_hcount) = kommundata(h_kommunindex(m_hcount)).kod
' h_la_region(m_hcount) = kommundata(h_kommunindex(m_hcount)).LA_region
h_BB_region(m_hcount) = BB_Region(h_kommunkod(m_hcount))
h_abroad(m_hcount) = h_abroad(oldhidx)
h_n_adults(m_hcount) = 1
h_size(m_hcount) = 1
h_n_child(m_hcount) = 0
h_emig_year(m_hcount) = i_emig_year(i)
h_wealth_debt_interest_varcomp_p(m_hcount) = 0 ' generates new random effects coefficient
h_wealth_debt_interest_varcomp_r(m_hcount) = 0 ' generates new random effects coefficient
h_wealth_debt_varcomp(m_hcount) = 0 ' generates new random effects coefficient
h_wealth_financial_varcomp(m_hcount) = 0 ' generates new random effects coefficient
h_wealth_debt_GLS1vc(m_hcount) = 0
h_wealth_debt_GLS2vc(m_hcount) = 0
h_wealth_debt_probit1vc(m_hcount) = 0
h_wealth_debt_probit2vc(m_hcount) = 0
'*** all housing-related variables are reset and recalculated in the housing module
h_house_owner(m_hcount) = 0
h_house_area(m_hcount) = 0
h_house_cost(m_hcount) = 0
h_house_interest(m_hcount) = 0
h_house_tax(m_hcount) = 0
h_wealth_real_home(m_hcount) = 0
h_wealth_real(m_hcount) = 0
h_new_housing(m_hcount) = 1 ' Marked for tenure choice
mark_h(m_hcount) = 0
'*** update information for old (female) household
h_n_adults(oldhidx) = 1
h_n_ge18(oldhidx) = h_n_ge18(oldhidx) - 1
h_size(oldhidx) = h_size(oldhidx) - 1
h_emig_year(oldhidx) = i_emig_year(indnr2index(inr_f))
h_wealth_debt_interest_varcomp_p(oldhidx) = 0 ' generates new random effects coefficient
h_wealth_debt_interest_varcomp_r(oldhidx) = 0 ' generates new random effects coefficient
h_wealth_debt_varcomp(oldhidx) = 0 ' generates new random effects coefficient
h_wealth_financial_varcomp(oldhidx) = 0 ' generates new random effects coefficient
h_wealth_debt_GLS1vc(oldhidx) = 0
h_wealth_debt_GLS2vc(oldhidx) = 0
h_wealth_debt_probit1vc(oldhidx) = 0
h_wealth_debt_probit2vc(oldhidx) = 0
'*** Handling of wealth variables according to specification by
'*** Anders Klevmarken 050320
' 1) Woman keeps residence (no action needed)
' 2) Woman takes over husbands postponed house taxes
i_housetax_postponed(indnr2index(inr_f)) = i_housetax_postponed(indnr2index(inr_f)) + _
i_housetax_postponed(i)
' 3) Woman takes over part of household debt
temp1 = h_wealth_debt(oldhidx)
h_wealth_debt(oldhidx) = maxi(0.5 * h_wealth_debt(oldhidx), _
mini(h_wealth_debt(oldhidx), 0.9 * h_wealth_real_home(oldhidx)))
h_wealth_debt(m_hcount) = temp1 - h_wealth_debt(oldhidx) ' Husbands share
' 4) Women keeps part of households financial wealth
' NOTE: adult children keeps part of financial wealth
' NOTE: if negative, the debt is increased accordingly
If h_n_ge18(oldhidx) - h_n_adults(oldhidx) > 0 Then
temp2 = 0.5 * h_wealth_financial(oldhidx) * 0.8 - 0.5 * h_wealth_real_home(oldhidx) + _
mini(temp1, 0.9 * h_wealth_real_home(oldhidx)) + 0.3 * i_housetax_postponed(i)
h_wealth_financial(m_hcount) = h_wealth_financial(oldhidx) - temp2
h_wealth_financial(oldhidx) = temp2
Else
temp2 = 0.5 * h_wealth_financial(oldhidx) - 0.5 * h_wealth_real_home(oldhidx) + _
mini(temp1, 0.9 * h_wealth_real_home(oldhidx)) + 0.3 * i_housetax_postponed(i)
h_wealth_financial(m_hcount) = h_wealth_financial(oldhidx) - temp2
h_wealth_financial(oldhidx) = temp2
End If
If h_wealth_financial(oldhidx) < 0 Then
h_wealth_debt(oldhidx) = h_wealth_debt(oldhidx) - h_wealth_financial(oldhidx)
h_wealth_financial(oldhidx) = 0
End If
If h_wealth_financial(m_hcount) < 0 Then
h_wealth_debt(m_hcount) = h_wealth_debt(m_hcount) - h_wealth_financial(m_hcount)
h_wealth_financial(m_hcount) = 0
End If
i_housetax_postponed(i) = 0
' 5) Other real wealth is divided 50/50
h_wealth_real_other(m_hcount) = h_wealth_real_other(oldhidx) / 2
h_wealth_real_other(oldhidx) = h_wealth_real_other(m_hcount)
' Allocate more memory for household index vectors if necessary
' NOTE: adding 100 elements at a time increases speed of computations
If UBound(hhnr2index) < största_hhnr Then
ReDim Preserve hhnr2index(1 To största_hhnr + 100) As Long
' ReDim Preserve hhnr2index_lag(1 To största_hhnr + 100) As Long
End If
hhnr2index(största_hhnr) = m_hcount
lifehist.write_hist i_indnr(i), "New hhnr (Leave women)"
End If
mark_i(i) = 0
Next
' keeping household vectors packed helps preventing errors
Call dyn_vect_h(m_hcount)
End Sub