Online Appendix for: Who Migrates and Why? Evidence from Italian Administrative Data Cristian Bartolucci

Claudia Villosio

Mathis Wagner

Collegio Carlo Alberto

Laboratorio R. Revelli

Boston College

December 2016

A

Return Migration

A considerable number of migrants return to their source region. We show that it is possible to analyze their return migration decision in the same framework as the original migration decision. Consider an individual i who has already migrated to the destination country and now every period t has the choice whether to stay Rit = 1 in the destination country j = n or return migrate Rit = 0 to the original source country j = s. As before, we assume she makes that decision based on the difference in outcomes yijt in each country and a one-time return migration cost rit . The return migration decision is carried out according to  Rit =

1 ”stay” iff yitn − yits − rit ≥ 0, 0 ”return migrate” iff yitn − yits − rit < 0.

As before, we will allow return migration costs rit to be a function of the time-varying x and time invariant α individual characteristics that affect outcomes (we omit the s superscripts), as well as other individual characteristics z that are excluded from the outcome equations, and unobserved individual-specific time-varying factors, which are 1

assumed to enter additively separable. Then we can rewrite the selection equation for potential return migrants as Rit = 1 (r (αi, xit , zit ) + ω it ≥ 0) ,

(1)

where ω it is distributed independently of α, x , and z with mean zero and variance σ 2ω . To incorporate the possibility of return migration into our basic model of migration we make the crucial assumption that the same potentially time-varying unobserved factors affect the return migration and migration decisions such that ω it ≡ vit . To provide a better idea for the intuition behind this restriction note that vit = uitn − uits − ζ M it , ω it = uitn − uits − ζ R it , where ζ M it are idiosyncratic factors that affect the decision to migrate from source to destination country, and ζ R it are idiosyncratic factors that affect the decision to stay in R the destination country and not return migrate. We assume that ζ M it ≡ ζ it , implying that

these factors are not moving costs per se, but rather affect the preference for living in a certain country. We can then combine selection equations for the migration decision, see Section 2 of the paper, and the return migration decision (see equation 1) into a selection equation describing whether the individual chooses to work in the destination country Nit = 1 or the source country Nit = 0: Nit = 1 (m (αi, xit , zit ) Di,t−1 + r (αi, xit , zit ) (1 − Di,t−1 ) + vit ≥ 0) ,

(2)

where Dit = 1 (yit∗ = yits ), i.e. if the individual was observed working in the source country last period. Then the outcome equations conditional on actually being observed contain two control

2

functions, one for migrants and one for return migrants: yitn |N =1 y its |N =0

  σN φ (m (αi , xit , zit )) φ (r (αi , xit , zit )) uv Dit + (1 − Dit ) +εitn , σv Φ (m (αi , xit , zit )) Φ (r (αi , xit , zit ))   σS φ (m (αi , xit , zit )) φ (r (αi , xit , zit )) = αi +β s xit − uv Dit + (1 − Dit ) +εits , σv 1 − Φ (m (αi , xit , zit )) 1 − Φ (r (αi , xit , zit ))

= ραi +β n xit +

where the coefficient on the control functions for potential migrants and potential return migrants are the same. Note that the set of observables can include variables that depend on whether you have been a migrant or not, for example, years spent in destination country. Hence, (return) migration can affect outcomes on account of differences in factor prices between the source and destination country, the returns to sorting, and because of differences in observables that are a function of the (return) migration decision. Intuitively, the key difference between the models with and without return migration is in the interpretation of the time-varying idiosyncratic factors. In the simpler model without return migration they are time-varying factors that affect the decision whether to migrate, essentially some type of moving cost. Once we allow for return migration in this framework they become time-varying location preferences. These might include personal, idiosyncratic events like the death of a parent, or a relationship break-up, as well as a broader socio-economic factors that differentially affect individuals.

B

Monte Carlo Studies

In our Monte Carlo study we present a simplified version of the model used in the empirical section. The model design is:

3

yit

  x β + α + u , if y ∗ > 0 it S i its =  x β + ρα + u , if y ∗ ≤ 0 it N



uits

   uitn  vit

i

itn

yit∗ = xit θ + ηzit + αi γ + vit    0 σ 0 σ Suv    S     ∼ N  0 ,  0 σN σN uv    S N 0 σ uv σ uv σ v

    

αi ∼ N (0, σ α ), cov(xit , αi ) 6= 0 N = 1000, T = (5, 10, 15, 20) γ 0 = 1, β S = 2, β N = 2, θ = 3, η = 5, ρ = 0.5, σ Suv = σ N uv = 0.5. We present results for estimates in 200 samples of 1000 individuals each. We describe the evolution of the estimators when T grows, reporting results for T = 5, 10, 15 and 20. Table 1 gives the Monte Carlo results for the estimators of γ, β S , β N , ρ, θ and η. In Figure 1, we show the performance of our estimation strategy in recovering the distribution of γ, β S , β N , ρ, θ and η. The strategy does well in recovering the parameters of the observable characteristics in the selection equation (θ and η), irrespective of the size of the longitudinal dimension. We find a significant difference in performance between the non bias-corrected and the bias-corrected estimators in recovering the coefficient on ability, γ, in the selection equation. The distribution of the bias-corrected estimate of γ is centered much closer to the true value than the non-corrected one. There is no significant improvement between T = 15 and T = 20, which is reassuring given the longitudinal dimension of the panel used in this paper. In Figure 1, we also present the performance of our estimation strategy in recovering the parameters that characterize the outcome equation. As expected, the coefficient on the observable characteristics perform well, while we find a difference in performance between the non bias corrected and the bias corrected coefficient of the unobserved fixed heterogeneity in the outcome equation in the North (ρ). The distribution of the bias corrected estimate of γ is closer to being centered around the true value, than the non

4

corrected one. Once again, there is no significant improvement between T = 15 and T = 20.

C

Full Information Maximum Likelihood Estimator

For the full information maximum likelihood estimator we consider the following model:

wit

  X β + α + ε , if s > 0 it S i it it =  X β + ρα + u , if s ≤ 0 it N

i

it

(3)

it

sit = Xit θ + Zit η + γαi + vit Together with the assumptions: 





σ 2ε

0 σ εv 0 0 ε   it      0  uit  0 σ 2u σ uv 0  ˜N  ,      0 σ εv σ uv σ 2v 0  vit     0 0 0 σ 2α 0 αi

       

This model can be estimated by maximum likelihood. The likelihood function is complicated but feasible. Define:

Θ = β S , β N , ρ, θ, η, γ, σ 2ε , σ 2u , σ 2v , σ uv , σ εv


0

The individual contribution to the likelihood is given by:   Z    α T Λi (Θ, ϕ) = log Πt=1 f (wit , sit |Xit , Zit , α, Θ) φ( )dα   σα R

and f (yit |xit , zit , α, Θ) is :

f (wit , sit |Xit , Zit , α, Θ) =

      

 φ

 w −X β −α  it it S i σε

wit −Xit β S −αi σ2 ε 2 1−(σ εv /σ 2 ε)

Xit θ+Zit ξz +γα+(σ εv /σ 2 ε)

 Φ

   if sit > 0

     w −X β −ρα    Xit θ+Zit η+γαi +(σ uv /σ 2  it it N i  u )wit −Xit β N −ραi  φ 1−Φ if sit ≤ 0 2 2 σ u

(1−σ uv /σ u )

5

Therefore the maximum likelihood estimators are:

  Z    α T N [   Π f (w , s |X , Z , α, Θ) φ( Σ log Θ = arg max )dα it it it it ML t=1 i=1 Θ   σα 

R

Note that FIML is computationally enormously more demanding than our proposed estimator. The contribution of each individual to the likelihood has an integral that has to be solved numerically. Since the integrand is normally distributed, the integral can be calculated numerically by Gaussian quadrature. The estimation involves calculating numerically N = 31, 626 integrals for each evaluation of the likelihood and each step in a standard Newton-Raphson Algorithm with 40 parameters requires 2 ∗ 40 + 2 ∗ 402 = 3280 evaluations. Moreover, this procedure has problems of convergence due to flat regions of the maximum likelihood function. Therefore, to avoid estimating derivatives numerically, we use the simplex algorithm for its optimization. In order to have a computationally feasible procedure we significantly reduce the number of observation keeping only individuals who are observed throughout the sample period. We recognize that the balanced panel contained in this sample is potentially selected and therefore, for comparison, we also re-estimate results using our iterative procedure. Results are presented in Table 2 using income as the outcome variable of interest. Most importantly, the FIML estimates are close to those of our iterative procedure in the balanced sample and, reassuringly, similar to those in the full sample. They also suggest negative selection on ability. The point estimate on ability is negative in the selection equation and the return to ability less than one in the selection equation.

D

Autocorrelated Shocks to Outcomes

Our fixed effect estimator with a selection correction is also consistent in the case of heteroscedasticity or serial correlation. We show Monte Carlo experiments of the model presented in Appendix B, above, adding serial dependence in the ui,t,s and ui,t,n . This serial correlation takes the following form: 6

ui,t,s ui,t,n

σ t,t−1 uit−1,s + ω i,t,s , =

σ t,t−1 ui,t−1,n + ω i,t,n ,

where ω i,t,s and ω i,t,n are N (0, 1) white noises. We estimate the coefficients in 200 samples of 500 individuals each observed for 10 periods. The primitives of the model are estimated with the iterative algorithm described in Section 4 of the paper and corrected with the jackknife bias correction. In figure 2, we present the distribution of the coefficients. We describe the evolution of the estimates when σ t,t−1 increases, reporting results for σ t,t−1 equal to 5%, 10%, 15% and 20%. We find that, as before, the coefficients are slightly biased due to the incidental parameter problem and that the jack-knife correction has a good performance reducing this bias. As expected, we do not observe major differences when we increase the degree of serial correlation. A second concern in the presence heteroscedasticity or serial correlation is the estimation of the asymptotic variance of the estimates. We obtain our standard errors and p-values by bootstrap, where sampling is over individuals, thereby dealing with serial correlation within individual.

E

Sample Selection

While the main sample used in the paper is representative of the full wage distribution in Southern Italy, Table 3 reports results for alternative samples: Panel A for wages and Panel B for income. In column one of the table, for convenience, we present again our baseline results. In column two, we limit the sample to those individuals who are observed working for at least five years in the South. This should decrease the estimation error in the ability term and thereby reduce the incidental parameter problem. In column three, we report results excluding all observations that correspond to return migration experiences. This is in the spirit of the extensive reduced-form literature, especially on 7

Mexico-US migration, which tends to focus on a single migration experience. Evidence for Mexico also finds different migrant selection patterns from urban and rural areas, see Fernandez-Huertas (2013). In columns four and five we present results separately for rural and urban Southern Italy, respectively. Limiting the sample to individuals who are observed in the South for at least five years reduces the number of observed migration events by 28 percent. Unsurprisingly, standard errors are larger but the results remain qualitatively the same. The degree of negative selection on ability is somewhat less pronounced than in the full sample, but remains statistically significant. This is reassuring, suggesting that our results are not driven by biases introduced by the incidental parameter problem. Dropping all observations associated with return migration experiences also does not result in any qualitative changes in our results. Our results do not seem to be substantially affected by the parsimonious way in which we model return migration decisions. Finally, the differences in selection between rural and urban areas are not pronounced.

8

Figure 1: Monte Carlo Simulations for Different Number of Observations per Individual (T ) bS with T=5

1.995

bS with T=10

2

2.005

1.996

1.998

bN with T=5

1.99

1.995

2

2

.48

2.005

.49

2.01

2.015

1.995

2

3

.5

.51

4

.52

.48

4

6

2

2.5

1

8

1.5

2.01

1.996

1.998

.5

3

10

4

.51

6

.485

.49

4

4.5

.6

.8

1

2

2.002

.495

2.5

2.002

1.998

1.999

2.006

1.996

1.998

8

4

.505

.51

.485

4

1.4

.6

.8

2.003

2.002

2.004

.495

.5

.505

.51

3

3.5

4

h with T= 20

6

1

2

.49

2.5

7

4

g with T=15

1.2

2.002

q with T= 20

3.5

5

2.001

r with st,t-1 = 20%

h with T=15

7

2

bN with T= 20

2.004

.5

3

g with T=10

2

2.001

q with T=15

3.5

5

2

r with T=15

h with T=10

g with T=5

.5

1.999

q with T=10

5

6

1.998

bS with T= 20

bN with T=15

2.005

.49

h with with T=5

2

2.004

r with T=10

q with T=5

2

2.002

bN with T=10

r with T=5

.47

bS with T=15

4.5

5

5.5

6

6.5

g with T= 20

1.2

1.4

.8

1

1.2

1.4

Note: Solid lines are the bias-corrected coefficients and dashed lines are the bias non-corrected ones. Monte Carlo experiments are based on 100 replications of a sample of 1,000 individuals for T = 5,10,15,20. Kernel = epanechnikov, bandwidth = 0:05. The vertical lines represent the true values of the parameters.

9

Figure 2: Monte Carlo Simulations for Different Degrees of Serial Correlation σ t,t−1 bS with st,t-1 = 5%

1.996

1.998

2

bS with st,t-1 = 10%

2.002

2.004

1.996

bN with st,t-1 = 5%

1.995

2

2.005

1.998

.49

.5

2.01

1.995

2

4

.51

6

.52

.48

.49

5

8

10

2

1

2

1.998

2.01

.51

6

15

0

5

1.995

.48

2

.49

2

3

.5

1

1.5

2

1.996

1.998

2.01

.51

6

0

5

1.995

.48

2

.5

1

1.5

2

2.01

.5

.51

.52

2

4

6

8

10

h with st,t-1 = 20%

15

2.5

2.005

.49

0

5

g with st,t-1 = 15%

3

2.004

q with st,t-1 = 20%

10

10

2.002

r with st,t-1 = 20%

.52

8

2

bN with st,t-1 = 20%

h with st,t-1 = 15%

15

2.5

2.005

.5

4

g with st,t-1 = 10%

2.5

2.004

q with st,t-1 = 15%

10

10

2.002

r with st,t-1 = 15%

.52

8

2

bS with st,t-1 = 20%

bN with st,t-1 = 15%

h with st,t-1 = 10%

10

1.5

2.005

.5

4

g with st,t-1 = 5%

.5

1.996

q with st,t-1 = 10%

h with st,t-1 = 5%

0

2.004

r with st,t-1 = 10%

q with st,t-1 = 5%

2

2.002

bN with st,t-1 = 10%

r with st,t-1 = 5%

.48

2

bS with st,t-1 = 15%

10

15

g with st,t-1 = 20%

3

.5

1

1.5

2

2.5

3

Note: Solid lines are the bias-corrected coefficients and dashed lines are the bias non-corrected ones. MonteCarlo experiment based on 100 replications of a sample of 500 individuals for σ t,t−1 = 5%, 10%, 15%, 20%. Kernel = epanechnikov, bandwidth = 0:05. The vertical lines represent the true values of the parameters.

10

Table 1: Monte Carlo Experiments

T =5

T = 10

T = 15

T = 20

γ βS βN ρ θ η γ βS βN ρ θ η γ βS βN ρ θ η γ βS βN ρ θ η

M ean .930 1.9994 2.0093 .4739 3.119 5.167 .935 1.9998 2.0055 .4835 2.968 4.932 .9588 2.0002 2.0049 .4859 2.979 4.955 .9628 1.9999 2.0046 .4872 2.973 4.946

Non-Corrected Coefficients M edian Std.Dev Q75 − Q25 .923 .174 .226 1.9992 .0029 .0036 2.0088 .0051 .008 .4736 .0109 .0137 3.032 .504 .617 5.029 .838 1.072 .923 .109 .157 1.9999 .0020 .0022 2.0059 .0032 .0049 .4829 .0075 .0099 2.915 .299 .369 4.844 .495 .615 .9557 .0889 .118 2.0003 .0014 .002 2.0048 .0024 .003 .4855 .0064 .0084 2.958 .237 .332 4.913 .400 .557 .9636 .0778 .103 1.9999 .0012 .0013 2.0045 .0022 .0028 .4875 .0051 .0062 2.961 .2216 .2958 4.920 .365 .4073

M ean .994 1.9994 2.0055 .4844 3.234 5.368 .980 1.9999 2.0028 .4916 3.047 5.075 .995 2.0002 2.0027 .4926 3.048 5.075 .9961 1.9998 2.0028 .4926 3.037 5.060

Bias-Corrected Coefficients M edian Std.Dev Q75 − Q25 .971 .234 .280 1.9992 .0030 .0035 2.0059 .0061 .0093 .4844 .0133 .0191 3.172 .677 .767 5.300 1.127 1.278 .953 .130 .179 1.9998 .0020 .0028 2.0032 .0040 .0049 .4909 .0092 .012 2.951 .359 .444 4.915 .599 .737 .998 .106 .146 2.0003 .0014 .0018 2.0025 .003 .0036 .4922 .0081 .0092 3.045 .277 .367 5.058 .469 .654 .9925 .0841 .1168 1.9999 .0012 .0013 2.0027 .0029 .0031 .4928 .0072 .0104 3.008 .234 .344 4.987 .389 .543

Note: 100 replications of a sample of 1,000 individuals for T = 5, 10, 15 and 20. The values of the parameters used in the DGP are: γ = 1, β S = 2, β N = 2, θ = 3, η = 5,ρ = 0.5. We report descriptive statistics for the sample of non-corrected coefficients in the first 4 columns and descriptive statistics for the sample of Jacknife-corrected coefficients in the last 4 columns.

11

Table 2: Full Information Maximum Likelihood (FIML) Estimates With Returnees Without Returnees FIML Iterative FIML Iterative Coeff. P-Value Coeff. P-Value Coeff. P-Value Coeff. P-Value Income in the South Part-Time -0.336 <0.001 -0.338 <0.001 -0.338 <0.001 -0.337 <0.001 Tenure 0.015 0.011 0.026 <0.001 0.01 0.034 0.026 <0.001 0.001 0.413 -0.001 <0.001 <0.001 0.427 -0.001 <0.001 Tenure Sq. Tenure Truncated <0.001 <0.001 0.061 0.002 <0.001 <0.001 0.059 0.002 Blue collar worker -0.082 <0.001 -0.074 <0.001 -0.086 <0.001 -0.073 <0.001 0.355 0.011 0.352 0.030 0.349 <0.001 0.352 0.030 Manager 0.015 0.054 0.018 <0.001 0.014 0.112 0.017 0.002 Pot. Experience -0.012 0.348 -0.001 <0.001 -0.007 0.146 -0.001 <0.001 Pot. Experience Sq. Multi-region Firm 0.057 <0.001 0.073 0.062 0.055 <0.001 0.063 0.090 Inverse Mills Ratio .261 <0.001 0.278 0.253 0.266 <0.001 0.274 0.156 Constant 0.407 <0.001 0.410 <0.001 0.403 <0.001 0.411 <0.001 Income in the North Part-Time -0.566 <0.001 -0.566 0.008 -0.645 <0.001 -0.641 0.010 Tenure 0.059 <0.001 0.064 <0.001 0.067 <0.001 0.065 0.004 Tenure Sq. 0.01 <0.001 -0.003 <0.001 0.003 <0.001 -0.003 0.002 Tenure Truncated -0.119 <0.001 -0.128 0.194 -0.169 <0.001 -0.192 0.118 Blue collar worker -0.284 <0.001 -0.283 0.006 -0.366 <0.001 -0.385 0.002 Manager -0.612 0.348 -0.612 0.162 -0.646 0.191 -0.338 0.106 Pot. Experience 0.002 0.489 0.003 0.393 -0.004 0.371 0.026 0.311 Pot. Experience Sq. 0.003 0.185 <0.001 0.335 -0.003 0.056 -0.001 0.403 Multi-region Firm 0.151 <0.001 0.148 0.168 0.256 0.011 0.061 0.096 Inverse Mills Ratio 0.142 <0.001 0.146 0.220 0.245 <0.001 -0.074 0.158 Ability (ρ) 0.352 <0.001 0.348 0.014 0.134 <0.001 0.352 0.072 0.059 0.059 0.054 0.379 -0.231 <0.001 0.018 0.397 Constant Selection Equation Part-Time -0.114 <0.001 -0.116 0.275 -0.215 <0.001 -0.217 0.206 -0.047 <0.001 -0.041 0.018 -0.044 <0.001 -0.033 0.080 Tenure Tenure Sq. -0.011 0.043 0.002 0.006 -0.010 0.011 0.002 0.028 Tenure Truncated -0.621 <0.001 -0.603 <0.001 -0.606 <0.001 -0.599 0.002 Blue collar worker -0.054 <0.001 -0.056 0.337 -0.119 <0.001 -0.120 0.236 Manager -1.23 <0.001 -1.23 <0.001 -1.26 <0.001 -1.26 <0.001 Pot. Experience -0.023 <0.001 -0.014 0.321 -0.014 <0.001 -0.012 0.792 Pot. Experience Sq. -0.005 0.011 -0.001 0.182 -0.001 0.112 -0.001 0.978 Multi-region Firm 0.87 <0.001 0.884 <0.001 0.981 <0.001 0.986 <0.001 Average Moves 0.453 <0.001 0.448 0.146 0.358 <0.001 0.370 0.132 No. Moves 0.094 <0.001 0.097 0.291 0.111 <0.001 0.121 0.405 Ability (γ) -0.084 <0.001 -0.083 0.373 -0.391 <0.001 -0.384 0.178 Contant 7.816 <0.001 7.829 <0.001 7.784 <0.001 7.793 <0.001 Note: These maximum likelihood estimates are for the model presented in equation (3). We only report results using income as the measure of labor market performance. Integral of α within each individual to the likelihood is estimated by a Gaussian Quadrature with 100 nodes. The likelihood function is optimized using the Simplex algorithm. Standard errors and p-values are obtained by bootstrap of the complete procedure with 100 iterations with replacement over individuals. The one-sided p-value is calculated as the proportion of sampled permutations where the sign of the estimated coefficient is different from the sign of the point estimate; see Cameron and Trivedi (2005, p. 363).

12

Table 3: Sample Selection (1) Baseline

Ability

Inverse Mills Ratio

Inverse Mills Ratio Ability

Ability

Inverse Mills Ratio

Inverse Mills Ratio Ability

(2) (3) (4) (5) At least 5 No Return Rural Urban years in South Migrants Panel A. Wages i. Selection Equation -0.257*** -0.138** -0.153** -0.302*** -0.223*** (0.041) (0.043) (0.059) (0.054) (0.041) ii. Wage Equation South 0.042** 0.027 0.042 0.057** 0.013 (0.021) (0.033) (0.024) (0.025) (0.018) iii. Wage Equation North 0.023* 0.044 0.121*** 0.01 0.031*** (0.030) (0.042) (0.101) (0.046) (0.033) 0.495*** 0.365*** 0.305*** 0.307*** 0.304 *** (0.030) (0.028) (0.036) (0.045) (0.041) Panel B. Income i. Selection Equation -0.362*** -0.214*** -0.334*** -0.317*** -0.419*** (0.017) (0.019) (0.021) (0.021) (0.024) ii. Income Equation South 0.984*** 0.104 0.942* 0.977*** 0.903*** (0.073) (0.093) (0.163) (0.113) (0.098) iii. Income Equation North 0.201*** 0.122** 1.203*** 0.136 0.265** (0.061) (0.093) (0.445) (0.087) (0.085) 0.441*** 0.081** 0.196 *** 0.117*** 0.153*** (0.026) (0.025) (0.072) (0.033) (0.040)

Note: Results are reported for the selection and outcome equations in the south and north of Italy. The reported coefficients are bias corrected, bootstrap standard errors are reported in brackets. * significant at 10%; ** significant at 5%; *** significant at 1%. Rural workers are defined as individuals start working in a Southern province with population density ≤ 200 inhabitants/KM2; urban workers are those start working in a Southern province with population density > 200 inhabitants/KM2. The ability term is the estimate of the fixed effect in the earnings (wages or income) equation in the South.

13