Kleiber - Applied econometrics in R

The output comprises the usual regression output along with the value of the log-likelihood and a Wald statistic paralleling the familiar regression F statistic. For convenience, a tabulation of censored and uncensored observations is also included. The results indicate that yearsmarried and rating are the main “risk factors”.

To further illustrate the arguments to tobit(), we refit the model by introducing additional censoring from the right:

R> aff_tob2 <- update(aff_tob, right = 4)

R> summary(aff_tob2)

Call:

tobit(formula = affairs ~ age + yearsmarried + religiousness + occupation + rating, right = 4, data = Affairs)

Scale: 7.94

Gaussian distribution

Number of Newton-Raphson Iterations: 4

Log-likelihood: -500 on 7 Df

Wald-statistic: 42.6 on 5 Df, p-value: 4.5e-08

The standard errors are now somewhat larger, reflecting the fact that heavier censoring leads to a loss of information. tobit() also permits, via the argument dist, alternative distributions of the latent variable, including the logistic and Weibull distributions.

Among the methods for objects of class “tobit”, we briefly consider a Wald-type test:

R> linear.hypothesis(aff_tob, c("age = 0", "occupation = 0"), + vcov = sandwich)

144 5 Models of Microeconometrics

Linear hypothesis test

Hypothesis: age = 0 occupation = 0

Model 1: affairs ~ age + yearsmarried + religiousness + occupation + rating

Model 2: restricted model

Note: Coefficient covariance matrix supplied.

Res.Df Df Chisq Pr(>Chisq)

Thus, the regressors age and occupation are jointly weakly significant. For illustration, we use a sandwich covariance estimate, although it should be borne in mind that, as in the binary and unlike the Poisson case, in this model, misspecification of the variance typically also means misspecification of the mean (see again Freedman 2006, for further discussion).

5.5 Extensions

The number of models used in microeconometrics has grown considerably over the last two decades. Due to space constraints, we can only a ord to briefly discuss a small selection. We consider a semiparametric version of the binary response model as well as multinomial and ordered logit models.

Table 5.2 provides a list of further relevant packages.

A semiparametric binary response model

Recall that the log-likelihood of the binary response model is

where F is the CDF of the logistic or the Gaussian distribution in the logit or probit case, respectively. The Klein and Spady (1993) approach estimates F via kernel methods, and thus it may be considered a semiparametric maximum likelihood estimator. In another terminology, it is a semiparametric singleindex model. We refer to Li and Racine (2007) for a recent exposition.

In R, the Klein and Spady estimator is available in the package np (Hayfield and Racine 2008), the package accompanying Li and Racine (2007). Since the required functions from that package currently do not accept factors as dependent variables, we preprocess the SwissLabor data via

R> SwissLabor$partnum <- as.numeric(SwissLabor$participation) - 1

which creates a dummy variable partnum within SwissLabor that codes nonparticipation and participation as 0 and 1, respectively. Fitting itself requires first computing a bandwidth object using the function npindexbw(), as in

R> library("np")

R> swiss_bw <- npindexbw(partnum ~ income + age + education +

+youngkids + oldkids + foreign + I(age^2), data = SwissLabor,

+method = "kleinspady", nmulti = 5)

A summary of the bandwidths is available via

R> summary(swiss_bw)

Single Index Model

Regression data (872 observations, 7 variable(s)):

Regression Type: Local-Constant

Bandwidth Selection Method: Klein and Spady

146 5 Models of Microeconometrics

Formula: partnum ~ income + age + education + youngkids + oldkids + foreign + I(age^2)

Objective Function Value: 0.5934 (achieved on multistart 3)

Continuous Kernel Type: Second-Order Gaussian

No. Continuous Explanatory Vars.: 1

Finally, the Klein and Spady estimate is given by passing the bandwidth object swiss_bw to npindex():

R> swiss_ks <- npindex(bws = swiss_bw, gradients = TRUE) R> summary(swiss_ks)

Single Index Model

Regression Data: 872 training points, in 7 variable(s)

Bandwidth: 0.383

Kernel Regression Estimator: Local-Constant

Confusion Matrix

Predicted

Actual 0 1

0 345 126

1 137 264

Overall Correct Classification Ratio: 0.6984 Correct Classification Ratio By Outcome:

0 1

0.7325 0.6584

McFadden-Puig-Kerschner performance measure from prediction-realization tables: 0.6528

Continuous Kernel Type: Second-Order Gaussian

No. Continuous Explanatory Vars.: 1

The resulting confusion matrix may be compared with the confusion matrix of the original probit model (see Section 5.2),

R> table(Actual = SwissLabor$participation, Predicted =

+round(predict(swiss_probit, type = "response")))

Predicted

Actual 0 1

no 337 134

yes 146 255

showing that the semiparametric model has slightly better (in-sample) performance.

When applying semiparametric procedures such as the Klein and Spady method, one should be aware that these are rather time-consuming (despite the optimized and compiled C code underlying the np package). In fact, the model above takes more time than all other examples together when compiling this book on the authors’ machines.

Multinomial responses

For illustrating the most basic version of the multinomial logit model, a model with only individual-specific covariates, we consider the BankWages data taken from Heij, de Boer, Franses, Kloek, and van Dijk (2004). It contains, for employees of a US bank, an ordered factor job with levels "custodial", "admin" (for administration), and "manage" (for management), to be modeled as a function of education (in years) and a factor minority indicating minority status. There also exists a factor gender, but since there are no women in the category "custodial", only a subset of the data corresponding to males is used for parametric modeling below.

To obtain a first overview of how job depends on education, a table of conditional proportions can be generated via

R> data("BankWages")

R> edcat <- factor(BankWages$education)

R> levels(edcat)[3:10] <- rep(c("14-15", "16-18", "19-21"),

+c(2, 3, 3))

R> tab <- xtabs(~ edcat + job, data = BankWages)

R> prop.table(tab, 1)

job

edcat custodial admin manage 8 0.245283 0.754717 0.000000 12 0.068421 0.926316 0.005263 14-15 0.008197 0.959016 0.032787 16-18 0.000000 0.367089 0.632911 19-21 0.000000 0.033333 0.966667

where education has been transformed into a categorical variable with some of the sparser levels merged. This table can also be visualized in a spine plot via

R> plot(job ~ edcat, data = BankWages, off = 0)

148 5 Models of Microeconometrics

Fig. 5.4. Relationship between job category and education.

or equivalently via spineplot(tab, off = 0). The result in Figure 5.4 indicates that the proportion of "custodial" employees quickly decreases with education and that, at higher levels of education, a larger proportion of individuals is employed in the management category.

Multinomial logit models permit us to quantify this observation. They can be fitted utilizing the multinom() function from the package nnet (for “neural networks”), a package from the VR bundle accompanying Venables and Ripley (2002). Note that the function is only superficially related to neural networks in that the algorithm employed is the same as that for single hidden-layer neural networks (as provided by nnet()).

The main arguments to multinom() are again formula and data, and thus a multinomial logit model is fitted via

R> library("nnet")

R> bank_mnl <- multinom(job ~ education + minority,

+data = BankWages, subset = gender == "male", trace = FALSE)

Instead of providing the full summary() of the fit, we just give the more compact

R> coeftest(bank_mnl)

z test of coefficients:

This confirms that the proportions of "admin" and "manage" job categories (as compared with the reference category, here "custodial") increase with education and decrease for minority. Both e ects seem to be stronger for the "manage" category.

We add that, in contrast to multinom(), the recent package mlogit (Croissant 2008) also fits multinomial logit models containing “choice-specific” (i.e., outcome-specific) attributes.

Ordinal responses

The dependent variable job in the preceding example can be considered an ordered response, with "custodial" < "admin" < "manage". This suggests that an ordered logit or probit regression may be worth exploring; here we consider the former. In the statistical literature, this is often called proportional odds logistic regression; hence the name polr() for the fitting function from the MASS package (which, despite its name, can also fit ordered probit models upon setting method="probit"). Here, this yields

R> library("MASS")

R> bank_polr <- polr(job ~ education + minority,

+ data = BankWages, subset = gender == "male", Hess = TRUE) R> coeftest(bank_polr)

z test of coefficients:

Estimate Std. Error z value Pr(>|z|)

using again the more concise output of coeftest() rather than summary(). The ordered logit model just estimates di erent intercepts for the di erent job categories but a common set of regression coe cients. The results are similar to those for the multinomial model, but the di erent education and minority e ects for the di erent job categories are, of course, lost. This appears to deteriorate the model fit as the AIC increases:

R> AIC(bank_mnl)

150 5 Models of Microeconometrics

[1] 249.5

R> AIC(bank_polr)

[1] 268.6

5.6 Exercises

1.For the SwissLabor data, plotting participation versus education (see Figure 5.1) suggests a nonlinear e ect of education. Fit a model utilizing education squared in addition to the regressors considered in Section 5.2. Does the new model result in an improvement?

2.The PSID1976 data originating from Mroz (1987) are used in many econometrics texts, including Greene (2003) and Wooldridge (2002). Following Greene (2003, p. 681):

(a)Fit a probit model for labor force participation using the regressors age, age squared, family income, education, and a factor indicating the presence of children. (The factor needs to be constructed from the available information.)

(b)Reestimate the model assuming that di erent equations apply to women with and without children.

(c)Perform a likelihood ratio test to check whether the more general model is really needed.

3.Analyze the DoctorVisits data, taken from Cameron and Trivedi (1998), using a Poisson regression for the number of visits. Is the Possion model satisfactory? If not, where are the problems and what could be done about them?

4.As mentioned above, the Affairs data are perhaps better analyzed utilizing models for count data rather than a tobit model as we did here. Explore a Poisson regression and some of its variants, and be sure to check whether the models accommodate the many zeros present in these data.

5.Using the PSID1976 data, run a tobit regression of hours worked on nonwife income (to be constructed from the available information), age, experience, experience squared, education, and the numbers of younger and older children.

152 6 Time Series

cross-sectional data. Dealing with time series data poses slightly di erent challenges. While we also need a rectangular, typically numeric, data matrix, in addition, some way of storing the associated time points of the series is required. R o ers several classes for holding such data. Here, we discuss the two most important (closely related) classes, “ts” and “zoo”.

R ships with the basic class “ts” for representing time series data; it is aimed at regular series, in particular at annual, quarterly, and monthly data. Objects of class “ts” are either a numeric vector (for univariate series) or a numeric matrix (for multivariate series) containing the data, along with a "tsp" attribute reflecting the time series properties. This is a vector of length three containing the start and end times (in time units) and the frequency. Time series objects of class “ts” can easily be created with the function ts() by supplying the data (a numeric vector or matrix), along with the arguments start, end, and frequency. Methods for standard generic functions such as plot(), lines(), str(), and summary() are provided as well as various time- series-specific methods, such as lag() or diff(). As an example, we load and plot the univariate time series UKNonDurables, containing the quarterly consumption of non-durables in the United Kingdom (taken from Franses 1998).

R> data("UKNonDurables")

R> plot(UKNonDurables)

The resulting time series plot is shown in the left panel of Figure 6.1. The time series properties

R> tsp(UKNonDurables)

reveal that this is a quarterly series starting in 1955(1) and ending in 1988(4). If the series of all time points is needed, it can be extracted via time(); e.g., time(UKNonDurables). Subsets can be chosen using the function window(); e.g.,

R> window(UKNonDurables, end = c(1956, 4))

Qtr1 Qtr2 Qtr3 Qtr4 1955 24030 25620 26209 27167 1956 24620 25972 26285 27659

Single observations can be extracted by setting start and end to the same value.

The “ts” class is well suited for annual, quarterly, and monthly time series. However, it has two drawbacks that make it di cult to use in some applications: (1) it can only deal with numeric time stamps (and not with more general date/time classes); (2) internal missing values cannot be omitted (because then the start/end/frequency triple is no longer su cient for