---
title: "Tools for Teaching Econometrics"
author: "Pierre Chausse^[University of Waterloo, pchausse@uwaterloo.ca]"
date: ""
bibliography: ref.bib 
output:
 pdf_document:
  number_sections: true
abstract: "This vignette explains how to use the different tools I used in my econometrics courses at the University of Waterloo. It includes functions to generate solutions from inference questions, to print regression results and more."
vignette: >
  %\VignetteIndexEntry{Tools for Teaching Econometrics}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteDepends{metricsUW}
  %\VignettePackage{metricsUW}
  %\VignetteEncoding{UTF-8}
header-includes:
- \DeclareMathOperator{\Ex}{E}
- \DeclareMathOperator{\var}{var}
- \DeclareMathOperator{\cov}{cov}
- \DeclareMathOperator{\pr}{Pr}
- \newcommand{\tp}[1]{#1^T}
- \newcommand{\reals}{\mathbb{R}}
---


# Introduction \label{sec:intro}

I use many R functions (@rbase) in my courses to print regression
results in a nice format, generate solutions for inference question on
the mean, the variance, and least squares models, to simulate data to
illustrate concepts, and so on. This document explains how to use them
using examples from the courses. To run the different functions
available in the package, you first need to load it:

```{r}
library(metricsUW)
```

# Printing regression results in equation format

To print a regression result from `lm` or `glm`, the method `printReg`
applied to `lm` and `glm` objects generates a latex equation with the
coefficient estimates and their standard errors. To have the equation
printed in Latex format in a R-Markdown document, the chunk option
`results='asis'` must be added. For example, the following is an
estimated wage equation using the `PSID1976` dataset from the `AER`
package (@aer):

```{r, results='asis'}
data(PSID1976, package="AER")
fit1 <- lm(wage~education+age+I(age^2), PSID1976)
printReg(fit1)
```

The function offers different options. You can add stars for
significant coefficients (`stars=TRUE`), adjusted $R^2$
(`adjrsq=TRUE`), limit the number of variables per line when the
equation does not fit the page (e.g. `maxpl=3`), replace default
standard errors by robust ones (`se=newse`), replace the default
t-distribution used to compute p-values by the N(0,1) distribution
(`dist="n"`) or omit variables. We present here a few examples.


- Adding adjusted $R^2$ and stars, and reducing the number of
  variables per line. I also increase the number of digits.

```{r, results='asis'}
fit2 <- lm(wage~education*city+heducation+age+I(age^2), PSID1976)
printReg(fit2, maxpl=3, adjrs=TRUE, stars=TRUE, digits=5)
```

- Replacing default standard errors by robust ones, and using the
  standard normal distribution for p-values (@sandwichlm).

```{r, results='asis'}
library(sandwich)
newse <- sqrt(diag(vcovHC(fit2)))
printReg(fit2, maxpl=3, se=newse, stars=TRUE)
```

For GLM estimation, the $R^2$ is replaced by residual deviance and AIC
is printed. Also, the left-hand side specifies that it is the link of
$\hat{Y}$ that has the linear representation. For example, the
following is the result from a Poisson regression with the log link,
using the dataset `fertil2` from the `wooldridge` package
(@wooldridge):

```{r, results='asis'}
data(fertil2, package="wooldridge")
fit3 <- glm(children~educ+age+I(age^2)+catholic+electric+radio+tv+heduc,
            family=poisson(link=log), data=fertil2)
printReg(fit3, maxpl=3, stars=TRUE)
```

If we just want to create a regression equation from a formula, the
method admits objects of class `formula`. Here is an example:

```{r, results='asis'}
printReg(log(wage)~education*female+age+I(age^2), maxpl=3)
```

# Solution to inference questions

The way the functions are organized is as follows. First the inference
function generates the solution, which is an object of class
`metricsSol`, and the `print` method generates the answer in Latex
format. It is meant to be printed inside an R-Markdown chunk, with
the option `results="asis"`. 

The purpose of these functions is to generate questions for
assignments, exercises or exams and have a fast way to produce a
detailed solution document. The questions can be generated by
datasets, in which case the solutions are derived from the same
dataset, or by made up numbers. 

## Introduction to Statistics

This section covers inference problems typically covered in an
introductory statistics course.

### Test on the mean

The solution is based on the following properties:

- If $X_i\sim N(\mu,\sigma^2)$, then $$
  \frac{\sqrt{n}(\bar{X}-\mu)}{\hat\sigma} \sim t_{n-1}\,,$$ where $n$
  is the sample size, $\bar{X}$ is the sample mean and $\hat\sigma$ is
  the sample standard errors defined as
  $\sum_{i=1}^n(X_i-\bar{X})^2/(n-1)$.
  
- If the distribution of $X_i$ is unknown, the true distribution of
  the test is also unknown, but the distribution it converges to is
  known. In my course, the t-distribution is only used when it is the
  exact distribution. If not, we use the N(0,1) as an
  approximation. The solution generator applies this rule as well.

Suppose we have a series and want to test the hypothesis $H_0:~\mu=c$
against $H_1:~\mu\neq c$, $H_1:~\mu>c$ or $H_1:~\mu<c$. The function
`testm` generates the solution. Consider the `wage` series form the
PSID1976 dataset, restricted to workers with positive hours:

```{r}
wage <- subset(PSID1976, hours>0)$wage
```

We want to test if the population average is equal to 5 dollars per
hour against the alternative that it is not equal to 5. We just insert
the following code in the chunk:

```{r, results='asis'}
testMean(wage, h0=5)
```

It is not necessary to use the `print` method directly, unless we want
to add something to the solution. For example, I like to add mark
distribution in my exam solutions:

```{r, results='asis'}
sol1 <- testMean(wage, h0=5)
print(sol1, addMess="1 point for the statistic and 1 point for the conclusion")
```

By default, we assume normality of the data, a size of 5\% and a
two-sided alternative. When we assume normality, the distribution used
for the critical value is the t-distribution with $(n-1)$ degrees of
freedom. If we don't, the critical value is based on the asymptotic
N(0,1) property of the test. It is also possible to change the
alternative hypothesis by either "greater" or "less" and the size of
the test:

```{r, results='asis'}
testMean(wage, h0=5, size=0.10, alter="less", assume="nonNormal")
```

By default, all decimals are kept to compute the statistic. This could
lead to slightly different solutions when the question is asked in an
exam, because the printed numbers are rounded. It is possible to round
the sample mean and standard errors before computing the test, through
the argument `dround`. For example, suppose an exam question was
generated directly from the data and the following table was printed
using `stargazer`:

```{r, results='asis'}
library(stargazer)
stargazer(data.frame(wage), digits=3, header=FALSE, float=FALSE)
```

We would generate the solution as follows:

```{r, results='asis'}
testMean(wage, h0=5, size=0.10, alter="less", dround=3)
```

This option exists for all solution generator, so we won't discuss it
further. It also possible to generate a solution without data. We just
need to provide the sample mean, the standard error and the sample
size:

```{r, results='asis'}
testMean(h0=4, xbar=3.7, se=0.8, n=40)
```

### Confidence interval for the mean

The $(1-\alpha)\times 100$\% confidence interval is defined as:

\[
[\bar{X} - q_{1-\alpha/2}\hat\sigma,~~\bar{X} + q_{1-\alpha/2}\hat\sigma]\,,
\]

were $q_{1-\alpha/2}$ is the $(1-\alpha/2)\times 100$\% quantile of
the $t_{n-1}$ when $X_i\sim N(\mu, \sigma^2)$, in which case the
coverage is exact, and the N(0,1) when the distribution of $X_i$ is
unknown. For the latter case, the coverage is just an approximation
based on the CLT. This rule inn consistent with the rule for tests on
the mean described in the previous section.

If we use the wage data from the previous section, the 95\% confidence
interval for the average wage, assuming normality, is:

```{r, results='asis'}
ciMean(wage)
```

As for `testMean`, we can choose not to assume normality and change
the `size`. It is also possible to round the mean and standard
deviation to match the solution of written questions.

```{r, results='asis'}
ciMean(wage, size=0.15, assume="nonNormal", dround=3)
```

Finally, we can specify the sample mean, standard deviation and sample size:

```{r, results='asis'}
ciMean(xbar=3.2, se=0.9, n=32)
```

### Tests on the difference between two means

The following properties are assumed in the solution generator
`testDiffMeans`. If we have two samples of sizes $n_1$ and $n_2$ for
$X_1$ and $X_2$, and want test the hypothesis $H_0:~\mu_1-\mu_2=c$, then:

- If $X_1\sim N(\mu_1, \sigma^2)$ and $X_2\sim N(\mu_2,\sigma^2)$, which implies that we assume equal variance, we have the following result under the null: $$\frac{\bar{X}_1-\bar{X}_2-c}{\hat\sigma\sqrt{\frac{1}{n_1}+\frac{1}{n_2}}}\sim t_{n_1+n_2-2}\,,$$ where $$\hat\sigma^2 = \frac{1}{n_1+n_2-2}\left[\hat\sigma^2_1(n_1-1) + \hat\sigma^2_2(n_2-1) \right]$$ and $\hat\sigma^2_1$ and $\hat\sigma^2_2$ are the usual bias corrected estimator of the variance of $X_1$ and $X_2$.

- For any other cases, which include non-normality and/or non-equal
  variances, the exact distribution of the test is unknown, so we use
  the approximated N(0,1) instead of the t-distribution. If the
  variances are not equal, the above test is not valid and must be
  replaced by:$$\frac{\bar{X}_1-\bar{X}_2-c}{\sqrt{\frac{\hat\sigma^2_1}{n_1} + \frac{\hat\sigma^2_2}{n_2}}}\approx N(0,1)\,.$$ 


> **Note**: The test presented here assumes that $Cov(X_1,X_2)=0$. We
> do not cover the non-zero covariance case in my courses, so it is
> not yet implemented.

Suppose we want to test if the average wage if the same for workers
living in a city and the ones not living in a city, we can proceed as
follows. We first consider the normal case with equal variances.

```{r, results='asis'}
wageCity <- subset(PSID1976, hours>0 & city=="yes")$wage
wageNoCity <- subset(PSID1976, hours>0 & city=="no")$wage
testDiffMeans(x1=wageCity, x2=wageNoCity, h0=0)
```

For any other cases, the approximated N(0,1) is used. Here is an
example with other specifications:

```{r, results='asis'}
testDiffMeans(x1=wageCity, x2=wageNoCity, h0=1, assumev="diff", size=0.10,
              alter="less")
```

As for the other tests, we can input estimated means and standard
errors instead of vectors. The arguments `xbar`, `se`, and `n` must be
vectors of two:

```{r, results='asis'}
testDiffMeans(h0=1, xbar=c(2.2,3.3), se=c(3.4,4.6), n=c(34,76))
```

### Confidence intervals for the difference between two means

The theory from the previous section also applies here: we use the
t-distribution only if the data is normally distributed and the
variances of $X_1$ and $X_2$ are the same. Also, the standard
deviation of $\bar{X}_1-\bar{X}_2$ is

\[
s= = \sqrt{\frac{1}{n_1+n_2-2}\left[\hat\sigma^2_1(n_1-1) + \hat\sigma^2_2(n_2-1) \right]
\left(\frac{1}{n_1}+\frac{1}{n_2}\right)}\,,
\]

if the variances are the same and the following if they are not:

\[
s = \sqrt{\frac{\hat\sigma^2_1}{n_1} + \frac{\hat\sigma^2_2}{n_2}}\,.
\]

The confidence interval for the difference in average wage between
workers from a city the ones not from a city, assuming normality and equal variance, is

```{r, results='asis'}
ciDiffMeans(wageCity, wageNoCity)
```

If we relax the normality and/or the equal variance we obtain:

```{r, results='asis'}
ciDiffMeans(wageCity, wageNoCity, assumev="diff", size=0.1)
```

As for testing the difference between means, we can also replace the
vectors of observations by `xbar`, `se`, and `n`.

### Test on the variance

The only implemented test at the moment is the one that assumes
normality. The assumption implies that under the null
$H_0:~\sigma^2=c$, we have the following distribution:

\[
\frac{(n-1)\hat\sigma^2}{c}\sim \chi^2_{n-1}
\]

Let's consider the following summary statistics from the `hprice1`
dataset:

```{r, results='asis'}
data(hprice1, package="wooldridge")
hprice1$lotsize <- hprice1$lotsize/1000
hprice1$sqrft <- hprice1$sqrft/1000
stargazer(hprice1[,1:5], digits=4, header=FALSE, float=FALSE)
```

Suppose we want to test $H_0:~\sigma^2=130$ against $H_1:~\sigma^2\neq
130$ for the lot size. In the following I show what happens if we set
`assume` to `"nonNormal"`. The test is performed, but a note is added
saying that the chi-square is not a valid distribution.

```{r, results='asis'}
testVar(hprice1$lotsize, h0=130, assume="nonNormal")
```

We can use one-sided tests and change the size:

```{r, results='asis'}
testVar(hprice1$lotsize, h0=130, alter="less", size=0.1)
testVar(hprice1$lotsize, h0=70, alter="greater", size=0.01)
```

We can also replace the numeric vector by values of `se` and `n`. For
example, we want to use the number from the table and test
$H_0:~\sigma^2=0.25$ for square footage, we proceed as follows:

```{r, results='asis'}
testVar(se=0.5772, n=88, h0=0.25)
```

### Confidence intervals for the variance

As for testing the variance in the previous section, the confidence
intervals presented here are only valid under normality. It is based
on the same chi-square distribution. For example, if we want to
generate a confidence interval for the variance of `sqrft` using the
data from the previous table, we can do it using the dataset (setting
`dround` to 4 as in the table):

```{r, results='asis'}
ciVar(x=hprice1$sqrft, dround=4)
```

This is equivalent to using the numbers from the table:

```{r, results='asis'}
ciVar(se=0.5772, n=88)
```

As for the test, we can change the `size` for a different confidence
level, and set `assume` to `"nonNormal"` to add a note about the
distribution not being valid under non-normality.

### Testing equality between two variances

We want to test $H_0:~\sigma^2_1=\sigma^2_2$ against
$H_1:~\sigma^2_1\neq \sigma^2_2$ at $\alpha\times 100$\%. The only
test I cover in the intro to statistics is the one that assumes
normality of the two samples and their independence. Under these
assumptions and the null hypothesis, we have the following
distribution:

\[
\frac{\hat\sigma^2_1}{\hat\sigma^2_2} \sim F(n_1-1, n_2-1)
\]

For this test to be two-sided, we define $\hat\sigma^2_1$ as the
highest estimate of the variance between the two samples. We reject of
the statistic exceeds the $(1-\alpha)\times 100$\% quantile of the
distribution. 

Consider the following summary statistics. The first table is for
female workers and the second is not male workers. The dataset used to
construct those tables is `wage1` form the `wooldridge` package.

```{r, results='asis', echo=FALSE}
data(wage1, package="wooldridge")
stargazer(subset(wage1, female==1)[,1:4], header=FALSE, float=FALSE)
stargazer(subset(wage1, female==0)[,1:4], header=FALSE, float=FALSE)
```

If we assume that wage is normally distributed and that female wage is
independent of male wage, we can test if the variance of wage is equal
between the two groups. First, using the dataset `wage1` with `dround=3` to
match the precision of the table:

```{r, results='asis'}
wF <- subset(wage1, female==1)$wage
wM <- subset(wage1, female==0)$wage
testDiffVar(wF, wM, dround=3)
```

We can also enter the information manually and change the size of the test:

```{r, results='asis'}
testDiffVar(se=c(2.529, 4.161), n=c(252, 274), size=0.1)
```

## Introduction to Econometrics

In this section, we present solution generators for inference based on
linear regressions. The general model is:

\[
Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + \cdots + \beta_kX_k + u
\]

By default, the estimated covariance matrix used to perform tests and
construct confidence intervals is the one obtained by `vcov`. When a
robust estimate is needed, it has to be obtained outside the solution
generator. As for inference on the means, the normality assumption
implies exact tests using either the t or the F distribution. Without
this assumption, the asymptotic standard normal or chi-square
distribution is used.

### Tests on a single regression coefficient

Consider the following estimated model:

```{r, results='asis'}
fit <- lm(wage~educ+exper, wage1)
printReg(fit)
```

If we want to test the hypothesis $H_0:~\beta_1=0.5$ against the
alternative $H_1:~\beta_1\neq 0.5$ with the assumption that the errors
are conditionally normal and homoskedastic, we can run the
following. Note that $\beta_1$ is the second coefficient, so we set
`which` to 2. Also, to generate a solution used the rounded
coefficient and standard error shown in the above regression, we set
`dround` to 4:


```{r, results='asis'}
testRegCoef(fit, which=2, h0=0.5, dround=4)
```

We can also obtain the solution without the estimated model, by
setting the coefficient estimate, its standard error, and sample size
and the number of coefficients to the values we see in the printed
equation:

```{r, results='asis'}
testRegCoef(h0=0.5, beta=0.6443, se=0.0538, n=526, ncoef=3)
```

If we relax the normality assumption, the critical value is obtained
using the N(0,1). If a vector of standard errors is provided, the
solution generator assumes that the ones obtained with `vcov` is
considered non-valid. In this case, there are no exact distribution,
so the N(0,1) is also used. In the following the test is performed
using the HCCM standard errors. For a solution based on a printed
regression, we have to print the regression with the robust standard
errors.

```{r, results='asis'}
rse <- sqrt(diag(vcovHC(fit)))
printReg(fit, se=rse)
```

In the following test, the N(0,1) is used instead of the t
distribution, because it is assume that the homoskedasticity
assumption does not hold.

```{r, results='asis'}
testRegCoef(fit, which=2, h0=0.5, dround=4, se=rse)
```

As for tests on the mean, we can change the size with the argument
`size` and the alternative with the argument `alter`.

### Confidence interval on a regression coefficient

The same rule applies for confidence intervals. It is exact (using the
t distribution) only if we assume normality and do not provide
alternative standard errors. For example, we can obtain the 90\%
confidence interval for the coefficient of `bdrms` in the following
regression as follows. First the estimated model is printed:

```{r, results='asis'}
fit <- lm(log(price)~bdrms+log(lotsize)+log(sqrft), hprice1)
printReg(fit)
```

Then we construct the interval:

```{r, results='asis'}
ciRegCoef(fit, which=2, size=0.10, dround=4)
```

To get confidence interval with robust standard errors or
non-normality, we just change `se` and `assume` as for the solution of
hypothesis tests.

### Testing linear combinations of coefficients (direct-t)

Any test on linear combinations of coefficients can be written as
$H_0:~R\beta=q$. In this section we consider single hypotheses, so $R$
has one row and $q$ is a scalar. For the function `testLinRegCoefs`, R
can be a plain vector or a matrix. The only requirement is to have a
length equal to the number of coefficients. Consider the following
regression, followed by the estimated covariance matrix:

```{r, results='asis'}
dat <- subset(wage1, female==1 & married==1 & smsa==1) 

fit <- lm(wage~educ+exper+tenure, dat)
printReg(fit, digits=5)
knitr::kable(round(vcov(fit),5))
```

If we want to generate a reproducible solution, we need to set
`dround` to 5, which is the number of digits shown in the regression
and covariance matrix. In the following, we test the hypothesis
$\beta_1-10\beta_2=0.6$, which implies $R=\{0,1,-10,0\}$ and $q=0.6$. 


```{r, results='asis'}
R <- c(0,1,-10,0)
q <- 0.6
testLinRegCoefs(object=fit, R=R, q=q, dround=5, alter="less")
```

If we set `assume` to `"nonNormal"` , or input an alternative
covariance matrix, the N(0,1) will be used to compute the critical
values.

```{r, results='asis'}
testLinRegCoefs(object=fit, R=R, q=q, dround=5, alter="greater",
                vcoef=vcovHC(fit))
```

It is also possible to generate a solution with specific covariance
matrix and coefficient vector. In this case, we need to specify if we
assume normality or not, in order to base the test on the t
distribution or the N(0,1). The following is the same solution without
the `lm` object, but with the t-distribution:

```{r, results='asis'}
v <- round(vcovHC(fit),5)
b <- round(coef(fit),5)
testLinRegCoefs(R=R, q=q, alter="greater", vcoef=v,  beta=b, n=nobs(fit))
```

### Testing multiple restrictions with the F test

The solution generator `testRegF` is for testing joint hypotheses
under the homoskedasticity assumption. I don't ask questions about the
non-homoskedastic case in exams, because it involves matrix
calculations that are not suite for questions to be answered in short
time. That's why this case is not covered. The test under normal errors is:

\[
\frac{(SSR_r-SSR_u)/q}{SSR_u/(n-k-1)} \sim F(q, n-k-1)\,,
\]

where $q$ is the number of restrictions being tested, $SSR_u$ is the
unrestricted SSR and $SSR_r$ is the restricted SSR. If we drop the
normality assumption, we use the following approximation based on the C.L.T.:

\[
\frac{(SSR_r-SSR_u)}{SSR_u/(n-k-1)} \approx \chi^2_q
\]

We do not present solutions based on the restricted and unrestricted
$R^2$, because its validity depends on left-hand side of the
regression being the same in the restricted and unrestricted
models. However, the solution of the joint significance of all
coefficients but the intercept is based on the $R^2$ of the
unrestricted model. We have the following with and without the
normality assumption:

\[
\frac{R^2/(k-1)}{(1-R^2)/(n-k-1)}\sim F(k-1, n-k-1)
\]
\[
\frac{R^2}{(1-R^2)/(n-k-1)}\approx \chi^2_{k-1}
\]

Consider the following models:

```{r, results='asis'}
data(hprice3, package="wooldridge")
hprice3 <- subset(hprice3, year==1978)
hprice3$price <- hprice3$price/1000
fitu <- lm(price~age+I(age^2)+rooms*baths+area+land, hprice3)
fitr <- lm(price~age+rooms+baths+area+land, hprice3)
printReg(fitu, maxpl=3, digits=5)
printReg(fitr, maxpl=3, digits=5)
```

We want to test the significance of all coefficients but the intercept
under normality. We can either provide the `lm` object of the
unrestricted model, or give the values of the $R^2$, sample size and
the number of coefficients:

```{r, results='asis'}
testRegF(fitu, dround=5)
testRegF(Rsq=0.62298, n=179, ncoef=8)
```

If we want to test the second model against the first, which
corresponds to $H_0:~\beta_2=\beta_7=0$, we proceed as follows
(`dround` is set to 2, because only two decimals of the SSR's are
printed). To see what happens if we do not assume normality, we set
`assume` to `"nonNormal"`.

```{r, results='asis'}
testRegF(fitu, fitr, dround=2, assume="nonNormal")
```

It is also possible to just provide the values of the SSR's and change the size:

```{r, results='asis'}
testRegF(assume="Normal", size=0.1, SSRu=62948.16, SSRr=66840.94,
         n=179, nrest=2, ncoef=8)
```

# Miscellaneous

In this section, I present some functions I used in the course as
tools for teaching or to generate random assignments. 

## Probability using a normal table

Let $X$ be distributed as $N(\mu,\sigma^2)$. In intro to statistics
courses, students have to compute $Pr(X\leq a)$, $Pr(X\geq a)$ or
$Pr(a\leq X\leq b)$ using a standardized normal table. The function
`probFromTable` explains using a latex equation, how to compute the
probability. The main arguments are `a`, `b` (if we want the last
probability), `type` (`"less"` or `"greater"` for the first two
probabilities), `mu` (E(X)), `sig2` (Var(X)), digits (the number of
digits shown in the normal table), `varName` (the name of the random
variable to be printed) and `approx` (is the normal an
approximation). By default, the function assumes that the normal table
shows $Pr(0\leq Z \leq z)$. Here is a few examples:

- $Pr(X\leq 8)$ for $X\sim N(10,2)$

```{r, results='asis'}
probFromTable(a=8, mu=10, sig2=2)
```

When generating a solution from exams or problem sets, we may be
interested in the numeric solution without the equation. This is done
by setting the argument `print` to `FALSE`:

```{r}
probFromTable(a=8, mu=10, sig2=2, print=FALSE)
```

- $Pr(X\geq 11)$ for $X\sim N(10,2)$

```{r, results='asis'}
probFromTable(a=11, type="greater", mu=10, sig2=2)
```

- $Pr(40 \leq Y \leq 48)$ for $Y\sim Binomial(100,0.5)$, which can be approximated by a $N(50,25)$. Note that `a` does not have to be the smallest number.

```{r, results='asis'}
probFromTable(a=48, b=40,  mu=50, sig2=25, approx=TRUE, varName="Y")
```

In some textbooks, the normal table shows the CDF instead of $Pr(0\leq
Z\leq z)$. If that's the case, we can set the argument `CDF` to
`TRUE`:

```{r, results='asis'}
probFromTable(a=48, b=40,  mu=50, sig2=25, approx=TRUE, varName="Y", CDF=TRUE)
```

## Generate statistics tables

For exercises or exam questions that require critical values from
statistics table, the function `distTable` created a Latex table using
the `xtable` package. For the t-distribution, it returns a table with
5 columns: 90\%, 95\%, 97.5\%. 99\% and 99.5\% quantiles. We can use
the default degrees of freedom, or enter them manually. By default, it
is a floating table with a title a the label `"tTable"` for
referencing. In the following, I set `float` to `FALSE`:

```{r, results='asis'}
distTable(df=c(4:9, Inf), dist="t", float=FALSE, digits=3)
```

For the chi-square table, the lower quantiles are also printed:

```{r, results='asis'}
distTable(df=c(4:10), dist="chisq", float=FALSE, digits=3)
```

The F-distribution has two sets of degrees of freedom, so a table only
contains one quantile at the time. We can change the default 95\% by
setting the argument size to something other than 0.05:

```{r, results='asis'}
distTable(df=c(4:10), df2=1:6, dist="f", size=0.10, float=FALSE, digits=3)
```

The last option is the Normal table. By default, it shows $Pr(0\leq
Z\leq z)$ for $Z\sim N(0,1)$. To have the CDF, we set `CDF` to `TRUE`.

```{r, results='asis'}
distTable(df=c(4:10), df2=1:6, dist="normal", CDF=TRUE,  float=FALSE, digits=3)
```

## Generate time series

The function `genTS` generates an `mts` object with two time series,
`X` and `Y`. One of the series is trend stationary (TS) and the second
one is a unit root process with possibly a constant or a linear
drift. The purpose of this dataset is to identify the process for each
series. The following is an example with all possible arguments, which
are all optional.

```{r, out.width="50%", fig.show="hold"}
dat <- genTS(n=200, frequency=4, start=1970, seed=112345)
plot(dat[,"X"], ylab="X")
plot(dat[,"Y"], ylab="Y")
```

Since we know the possible processes, it is not too hard to determine
the process of each series:

```{r, out.width="50%", fig.show="hold"}
fitX <- lm(X~time(X) + I(time(X)^2), dat)
Xdt <- ts(residuals(fitX), start=start(dat), frequency=frequency(dat))
fitY <- lm(Y~time(Y) + I(time(Y)^2), dat)
Ydt <- ts(residuals(fitY), start=start(dat), frequency=frequency(dat))
acf(Xdt, main="ACF of the detrended X")
acf(Ydt, main="ACF of the detrended Y")
```

We would conclude that since the detrended $Y$ is highly persistent,
it is probably a unit root process, and $X$ would be a trend
stationary series. The true processes are saved in the attribute
`"which"`:

```{r}
attr(dat, "which")
```

## Generate supply and demand functions

This tool is to illustrate the identification problem when estimating
demand and supply equations. The main function is `genDemSup`, which
generates a dataset that contains price, quantity, exogenous variables
and errors. The system is:

\begin{eqnarray*}
Q^d &=& d_1 + d_2P + d_3X + \eta\\
Q^s &=& a_1 + a_2P + a_3Z + e\\
\end{eqnarray*}


The function has the following arguments:

- `n`: The sample size
- `a`: The vector of supply coefficients.
- `d`: The vector of demand coefficients.
- `sige`, `sigeta`, `sigX`, `sigZ`: The variances of $e$, $\eta$, $X$ and $Z$.

It returns an object of class `DemSup`, with a `print` method that
generates the equations in Latex format, and a `plot` method that
plots the data points and a few demand and supply curves with their
equilibrium points. The coefficients must be such that an equilibrium
exists, which is the case of the default values. The following shows
the system of equation with by default:

```{r, results='asis', fig.align="center", out.width="50%"}
set.seed(1223345)
sys <- genDemSup(200)
print(sys)
```

The element `dat` of the object is a `data.frame` that contains $Q$,
$P$, $X$, $Z$, $e$ and $\eta$ (as `eta`). We can therefore try to estimate the demand function by OLS:

```{r, results='asis'}
fit <- lm(Q~P+X, sys$dat)
printReg(fit)
```

We could then plot some demand and supply with the fitted line to
illustrate the identification problem. The arguments of `plot` are

- `x`, the `DemSup` object
- `nCurves`: The number of demand and supply curves to plot. The
  default is all curves.
- `nPoints`: The number of data points to plot. The default is all
  points.
- `fac`: This is a vector of 2. It determines the range for each
  demand and supply curve plotted. The price range is
  $P^*\times$`fac[1]` to $P^*\times$`fac[2]`, where $P^*$ is the
  equilibrium price. By default it goes from $0.9P^*$ to $1.1P^*$.

In the following, we plot the fitted line when $X=\bar{X}$:

```{r, results='asis', fig.align="center", out.width="50%"}
plot(sys, nCurves=10)
pr <- predict(fit, newdata=data.frame(P=50:70, X=mean(sys$dat$X)))
lines(50:70, pr, col=2, lty=2, lwd=2)
```

We can use the data to analyze the exogeneity condition, since we are
given the vector of true errors. For example, if we want to estimate
the demand, we need to show that $Z$ is exogenous. It can be done by
regressing $Z$ on both errors and check the significance of the
f-test:

```{r, results='asis'}
fit2 <- lm(Z~e+eta, sys$dat)
printReg(fit2)
testRegF(fit2)
```

We therefore, conclude that $Z$ is exogenous, because it cannot be
explained by the error terms. We can also look at the relevance of the
instrument by estimating the reduced form equation and test the
significance of the coefficient of $Z$:

```{r, results='asis'}
fit2 <- lm(P~Z+X, sys$dat)
printReg(fit2, stars=TRUE)
```

Since the coefficient of $Z$ is highly significant, we conclude that
the instrument is relevant. We can then estimate the slope by
two-stage least squares and plot the fitted line:

```{r, results='asis', fig.align="center", out.width="50%"}
sys$dat$Phat <- fitted(fit2)
fit3 <- lm(Q~Phat+X, sys$dat)
printReg(fit3)
```

```{r, results='asis', fig.align="center", out.width="50%"}
plot(sys, nCurves=10)
lines(50:70, pr, col=2, lty=2, lwd=2)
pr2 <- predict(fit3, newdata=data.frame(Phat=50:70, X=mean(sys$dat$X)))
lines(50:70, pr2, col="darkgreen", lty=2, lwd=2)
legend("topleft", c("OLS", "TSLS"), col=c(2, "darkgreen"), lty=2, lwd=2, bty='n',
       cex=.7)
```

The parameter values can be use to illustrate different aspect of the
simultaneity problem. For example, we can show that the demand is
identified only the supply shifts (`sigeta=0` and `sigX=0')

```{r, results='asis', fig.align="center", out.width="50%"}
set.seed(1122345)
sys <- genDemSup(200, sigeta=0, sigX=0)
plot(sys, nCurves=10)
```

We can also illustrate the weak instrument problems. The `tsls`
function is from the `momentfit` package (@momentfit) and the `texreg`
function from the `texreg` package (@texreg). The `extract` method
needed to run `texreg` on `tsls` objects is included in the
`metricsUW` package.

```{r, message=FALSE, results='asis'}
library(texreg)
library(momentfit)
sysS <- genDemSup(200, a=c(20,1,-2))
sysW <- genDemSup(200, a=c(20,1,0.01))
fitS <- tsls(Q~P+X, ~X+Z, data=sysS$dat)
fitW <- tsls(Q~P+X, ~X+Z, data=sysW$dat)
texreg(list(Strong=fitS, Weak=fitW), digits=4, table=FALSE)
```

We see that the weakness of the instrument (represented here by
$a_3=0.01$ instead of $a_3=-2$) increase the variance of the
coefficients. The weakness is also clear from the first stage
F-test. Similarly, we can see what happens if the instrument has a
small variance by changing the value of `sigZ`. 

The tool also offers the possibility of having an over-identified
model. By default, we have `d=c(100, -1, 2, 0)` and `sigX=c(5,5)`,
which implies that only one exogenous variable is included in the
demand function. By replacing `d[4]` by a non-zero value, the demand
becomes an equation with the exogenous variables $X_1$ and $X_2$, and
`sigX` contains the variance of the two variables:

```{r, results='asis', out.width="50%", fig.align="center"}
sys2 <- genDemSup(n=200, d=c(100, -1, 2, 0.5))
print(sys2)
plot(sys2, nCurve=15)
```

We can test the strength of the instruments using the `testRegF` function:

```{r, results='asis'}
fitu <- lm(P~X1+X2+Z, sys2$dat)
fitr <- lm(P~Z, sys2$dat)
printReg(fitu)
printReg(fitr)
testRegF(fitu, fitr)
```

We therefore have strong instruments. We can also test the exogeneity
of the instruments, by using the `specTest` method of `momentfit` or
the `printReg` method for objects of class `tsls`:

```{r, results='asis'}
fit <- tsls(Q~P+Z, ~X1+X2+Z, data=sys2$dat)
printReg(fit, strength=TRUE, Jtest=TRUE, Exo=TRUE)
```

```{r}
specTest(fit)
```

# References