Exercices.Rmd

---
title: "Exercices, ECN 6338, Hiver 2023"
author: "William McCausland"
date: "`r Sys.Date()`"
output: pdf_document
urlcolor: blue
---

## Cours 1. Introduction

### Exercices théorique

1. Montrez que $6n^2 + 3n + 10 = O(n^2)$.

### Exercices de computation

1. Lisez l'aide sur la fonction \texttt{expm1} et démontrez son avantage par rapport à la fonction \texttt{exp} pour évaluer $e^x - 1$ quand $|x|$ est prés de zéro.
Vous pouvez suivre l'exemple sur \texttt{log1p} dans les diapos.

1. Téléchargez le paquet R \texttt{microbenchmark}, lisez l'aide sur la fonction \texttt{microbenchmark} et mesurer le temps nécessaire pour faire les opérations \texttt{x * y}, \texttt{x / y}, \texttt{x\textasciicircum y}, \texttt{exp(y)} et \texttt{log(x)}, pour un vecteur $x > 0$ de mille éléments et un vecteur $y$ de mille éléments.

\pagebreak

## Cours 2. La résolution de systèmes d'équations linéaires

### Exercices préliminaires

1. Soit $L_1$ et $L_2$ des matrices triangulaires inférieures $n\times n$.
Soit $U_1$ et $U_2$ des matrices triangulaires supérieures $n\times n$.
Supposez que $L_1$ et $U_1$ sont inversibles.
    a. Parmi les matrices suivantes, lequelles sont toujours triangulaires inférieures :
    $L_1L_2$, $L_1+L_2$, $L_1^{-1}$, $L_1^\top$, $U_1^\top$, $L_1U_1$?
    a. Parmi les matrices suivantes, lequelles sont toujours triangulaires supérieures :
    $U_1U_2$, $U_1+U_2$, $U_1^{-1}$, $U_1^\top$, $L_1^\top$, $L_1U_1$?
    
1. (Substitution avant et substitution arrière)
Soit $L$ et $U$ des matrices inversibles $n \times n$, où $L$ est triangulaire inférieure et $U$ est triangulaire supérieure.
Soit $y$ un vecteur $n \times 1$.
Notez que Judd (Algorithme 3.1, page 58) utilise le terme "back substitution" pour décrire deux algorithmes distincts (a. et b. ci-bas).
Plusieurs auteurs font une distinction entre "back substitution" (b.) et "forward substitution" (a.).
    a. Trouver un algorithme pour résoudre l'équation $Lx = y$.
    Notez que $L_{11} x_1 = y_1$, alors $x_1 = y_1/L_{11}$.
    a. Trouvez un algorithme pour résoudre l'équation $Ux = y$.
    Commencez par $x_n$.

1. Soit $x$ un $n$-vecteur aléatoire avec moyenne $\mu$ ($n\times 1$) et variance $\Sigma$ ($n \times n$). Soit $A$ une matrice constante $m \times n$.
Montrez que la moyenne et la variance de $Ax$ sont $A\mu$ et $A\Sigma A^\top$, pour le cas où $\Sigma$ est diagonal.

### Exercices théoriques

1. Quelquefois, il est plus facile de spécifier la matrice de précision $H = \Sigma^{-1}$ que la matrice de variance $\Sigma$, pour une loi gaussienne multivariée $N(\mu,\Sigma)$.
Si vous avez la matrice $H$ et non la matrice $\Sigma$, décrivez comment on peut tirer des variables aléatoires $N(\mu,\Sigma)$ et évaluer la densité gaussienne multivariée.
Commencez par la décomposition cholesky de $H$.

1. La matrice de précision $H$ pour un processus AR(1) gaussien est
tridiagonale ($H_{ij} = 0$ pour $|i-j|>1$) et symmétrique. Si $\omega = \sigma^{-2}$ est la précision de l'innovation et $\phi$ est le coefficient d'autorégression, le diagonal de $H$ est $\omega(1,1+\phi^2,\ldots,1+\phi^2,1)$ et le sous-diagonal est $-\omega(\phi, \ldots, \phi)$.
Trouvez la décomposition cholesky de $H$, en utilisant une description de l'algorithme (page 59 de Judd, par exemple).

1. Dans l'exemple MCO, calculez le vecteur $e = y - Xb$ des résiduelles et $\hat{\sigma}^2 (X^\top X)^{-1}$, une estimation de la variance de $b$. Utilisez seulement $y$ et la décomposition QR de $X$, pas $X$ elle-même. Notez que $\hat{\sigma}^2 = e^\top e/(n-k)$.

1. Donnez la matrice creuse des diapos en format colonne compressée après les deux opérations suivantes. Décrivez brièvement les algorithmes pour l'insertion et la supression en général.
    a. Une insertion \texttt{A[2,4] <- 24},
    a. Une supression \texttt{A[3,2] <- 0}.

### Exercices de computation

1. Regardez la démonstration \texttt{GPA.R} (au site Github du cours) de l'analyse d'une régression linéaire, exemple 3.1 dans Wooldridge (2016) Introductory Econometrics.
Calculez $b$, $e$ et $\hat{\sigma}^2 (X^\top X)^{-1}$ en utilisant la décomposition QR.
La fonction \texttt{qr} effectue la décomposition QR, les fonctions \texttt{qr.Q} et \texttt{qr.R} extraient les matrices $Q_1$ et $R_1$ du résultat.
Les fonctions \texttt{backsolve} et \texttt{t} (transpose) pourraient être utiles.

\pagebreak

## Cours 3. Quelques sujets préalables

### Exercices préliminaires

1. Trouvez le panier de consommation $(x_1^*,x_2^*)$ qui maximise $U(x_1,x_2) = x_1^{1/2} x_2^{1/2}$ sous les contraintes $p_1x_1 + p_2x_2 = m$, $x_1 \geq 0$ et $x_2 \geq 0$, où $m > 0$, $p_1 > 0$ et $p_2 > 0$.
Utilisez la méthode de Lagrange.

1. Si $X_i|\lambda \sim \mathrm{iid}\;\mathrm{Exp}(\lambda)$ (loi exponentielle) et $\lambda \sim \mathrm{Ga}(\alpha, \beta)$ (loi gamma)
écrivez
    a. la fonction de densité conditionnelle de $(X_1,\ldots,X_n)$ sachant $\lambda$.
    a. la densité conjointe de $(X_1,\ldots,X_n)$ et $\lambda$.

### Exercices théoriques

1. Considérez l'exercice préliminaire 1, avec la contrainte remplacée par $(p_1x_1 + p_2x_2 - m)^3 = 0$.
La solution du problème doit être pareille.
Décrivez le problème qui se présente quand vous essayez d'utiliser la méthode de Lagrange.
Expliquez brièvement pourquoi ce problème ne représente pas une contradiction du théorème de Lagrange.

1. Pour le modèle de l'exercice préliminaire 2,
    a. trouvez l'estimateur maximum de vraisemblance $\hat\lambda$ de $\lambda$; trouvez la moyenne et la variance de $\hat\lambda^{-1}$. (La moyenne et la variance de $\hat\lambda$ n'ont pas d'expression simple.)
    a. trouvez la distribution postérieure $\lambda|x^\circ$, où $x^\circ$ est une réalisation de $(X_1,\ldots,X_n)$.

### Exercices de computation

Aucune, étant donnée la nature du cours 3.

\pagebreak

## Cours 4. L'optimisation statique

### Exercice préliminaire

1. Lisez la section "Newton's Method" pp 96-97 dans Judd.
1. Prenez la fonction de log vraisemblance (graphique, diapo 19, cours 3) $L(\theta;y) = n_1 \log \theta + n_0 \log (1-\theta)$, pour $n_0 = 200$ et $n_1 = 230$.
    a. Calculez l'expansion Taylor d'ordre deux (quadratique) de $L(\theta;y)$ autour de $\theta = 1/2$.
    a. Trouvez la valeur de $\theta$ qui maximise cette fonction quadratique.
    a. Considerez le même problème que a., mais pour la vraisemblance ${\cal L}(\theta;y) = \theta^{n_1} (1-\theta)^{n_0}$.
    Faites assez de ce problème pour vous convaincre que c'est plus difficile.

### Exercices théoriques

1. On veut modeler le choix binaire $y_i \in \{0,1\}$ d'une individu $i$ comme fonction du vecteur $x_i \in \mathbb{R}^k$ des variables explicatives, pour un échantillon aléatoire $\{(y_i, x_i)\}$, $i=1,\ldots,n$.
Le modèle probit donne la probabilité conditionnelle $\Pr[y_i|x_i]$ comme $\Pr[y_i|x_i] = \Phi(x_i^\top \beta)$ où $\beta \in \mathbb{R}^k$ et $\Phi$ est la fonction de répartition d'une variable aléatoire gaussienne standard (\texttt{pnorm} en R; sa dérivée de première ordre est la densité, \texttt{dnorm} en R).
    a. Écrivez la fonction de log vraisemblance $L(\beta;y, X)$,
    où $y = (y_1,\ldots,y_n)$ est $X$ et la matrice $n \times k$ de variables explicatives.
    a. Donnez le gradient et la matrice hessienne de $L(\beta;y, X)$.
    Montrez que la matrice hessienne est négative définie.

### Exercices de computation

L'exercice est de reproduire les résultats du modèle probit de l'exemple 17.1 de Wooldridge, qui concerne le choix des femmes mariées d'être dans la population active ou non.
Vous pouvez charger les données avec
```{r mroz}
library(wooldridge)
data(mroz, package='wooldridge')
```
et ensuite voir les statistiques sommaires avec \texttt{summary(mroz)}.

1. Écrivez des fonctions en R pour calculer la valeur, le gradient et la matrice hessienne de la fonction de log vraisemblance $L(\beta;y, X)$ du modèle probit.
1. Trouvez l'estimation maximum de vraisemblance $\hat\beta$ pour le modèle probit de l'exemple 17.1 de Wooldridge, avec les méthodes Nelder-Mead et BFGS.
Les variables sont une constante,
\texttt{nwifeinc}, \texttt{educ}, \texttt{exper}, \texttt{exper} carré, \texttt{age}, \texttt{kidslt6}, \texttt{kidsge6}.
1. Calculer moins la matrice hessienne de la log vraisemblance au point $\hat\beta$.
Calculer l'inverse de cette matrice et donnez les écarts-types asymptotiques des éléments de $\hat\beta$ (les racines carrées des éléments diagonaux).



\pagebreak

## Cours 5. La résolution de systèmes d'équations non-linéaires

### Exercices préliminaires

1. Considérez le problème de trouver une racine de la fonction $f(x) = e^x - x - 2$ dans l'intervalle $[0,2]$. Faites avec une calculatrice :
    a. trois cycles de l'algorithme 5.1 (bisection).
    a. trois cycles de l'algorithme 5.2 (Newton).

### Exercices théoriques

1. Si on utilise la méthode de dichotomie pour trouver une racine de la fonction $f(x)$ dans l'intervalle $[0,1]$, combien de fois est-ce qu'il faut évaluer $f$ pour obtenir un intervalle de longeur $\sqrt{\epsilon}$ (ou moins) qui contient une racine, où $\epsilon$ est l'epsilon de la machine utilisée, égale à $2.22 \times 10^{-16}$.

1. Considérez une économie avec deux agents et deux biens.
Agent $a$ a une dotation de $(\tfrac{3}{4},\tfrac{1}{4})$ des deux biens et une
fonction d'utilité $u^a(x_1^a, x_2^a)$, où $(x_1^a, x_2^a)$ représente sa consommation des deux biens.
Agent $b$ a une dotation de $(\tfrac{1}{4},\tfrac{3}{4})$ et une fonction d'utilité $u^b(x_1^b, x_2^b)$.
L'agent $b$ peut offrir un échange à l'agent $a$, à prendre ou à laisser.
Supposons que $a$ accepte n'importe quelle échange où $a$ consomme un panier avec utilité égale ou supérieur à sa dotation (nonobstant les résultats empiriques sur les jeux de l'ultimatum).
Supposons que les deux fonctions d'utilité sont strictement croissantes et concaves, ainsi que continument dérivables.
Donnez un système $f(x)=0$ de deux équations en deux variables dont on peut utiliser la solution pour trouver le choix optimal de $b$.

### Exercices numériques

1. Pour le jeu de duopole (pages 162-163 de Judd et plusieurs diapos),
    a. écrivez des fonctions en R pour évaluer les fonctions $u(e^y,e^z)$, $u_Y(e^y,e^z)$, $u_Z(e^y,e^z)$,
    $\Pi^Y_1(e^y,e^z)$ et $\Pi^Z_2(e^y,e^z)$ et testez ces fonctions en utilisant
    les dérivées numériques.
    a. trouvez l'équilibre en utilisant la fonction \texttt{broyden} du paquet \texttt{pracma}.
    a. évaluez les gains et les pertes (pour les deux firmes et le consommateur) si firme 2 trouve une technologie qui réduit le coût de production unitaire de $Z$ de 0.08 à 0.06. 



\pagebreak

## Cours 6. L'approximation de fonctions

### Exercices préliminaires

1. Montrez que les trois premiers polynômes de Legendre, $1$, $x$ et $(3x^2-1)/2$
sont orthogonaux selon le produit intérieur défini à la page 203 de Judd.
Regardez le tableau 6.3 pour la spécification de $a$, $b$ et $w(x)$ pour les
polynômes de Legendre.

### Exercices théoriques

1. Trouvez l'approximation Padé (2,1) de $f(x) = \log x$ au point $x_0 = 1$.

1. Trouvez les valeurs des paramètres $\alpha$ et $\beta$ de la fonction de
pondération $w(x)$ associée aux polynômes Jacobi, pour les cas spéciaux
    a. Legendre
    a. Tchebyshev
    
1. Soit $f(x) = x^3$. Trouvez le polynôme $p_2$ d'ordre 2 qui minimise
$\langle f - p_2, f - p_2 \rangle$, pour le produit intérieur défini par
\[
  \langle f, g \rangle = \int_0^\infty e^{-x} f(x) g(x)\, dx.
\]
Notez que
\[
  \int_0^\infty x^k e^{-x} = k!.
\]

1. Montrez que $b_{i,n}'(x) = n(b_{i-1,n-1}(x) - b_{i,n-1}(x))$, où $b_{i,n}(x)$
est le $i$-ième polynôme Bernstein d'ordre $n$ et $b_{i,n}'(x)$ est sa dérivé.

### Exercices numériques

1. Considérez l'approximation de la fonction $f(x) = \sin(\pi x/2)$ sur
l'intervalle $[0,4]$, un cycle complet.
Faites les graphiques de la fonction et des approximations suivantes :
    a. Approximation Taylor d'ordre $n=3$ autour de $x_0 = 2$,
    a. Approximation Bernstein d'ordre $n=10$,
    a. Approximation spline cubique d'hermite qui utilise la fonction $f$ et sa
    première dérivé $f'$, évalués aux points $x=0,1,2,3,4$.

\pagebreak

## Cours 7. L'intégration et la dérivation (méthodes déterministes)

### Exercices préliminaires

1. Évaluez l'intégrale $\int_0^{\pi/2} \sin x\, dx$.
Lisez la page 260 de Judd et approximez l'intégral par quadrature
Gauss-Legendre avec $N=3$ points.

### Exercices théoriques

1. Supposez que $f \in C^4$.
Utilisez le théorème de Taylor-Lagrange pour trouvez la deuxième dérivée $f''(x)$ comme la somme des deux termes suivants :
    a. une formule, une fonction linéaire de $f(x)$, $f(x+h)$ et $f(x-h)$,
    a. un terme résiduel proportionnel à $h^2$.

### Exercices numériques

1. Dans l'exemple de la section 7.6 de Judd, trouvez la valeur optimale de $\omega$ pour les valeurs suivantes des paramètres : $\gamma = -5$, $\mu=0,2$, $\sigma=0,3$, $R=1,1$  $W=2$.
Utilisez l'option BFGS de la commande \texttt{optim} pour faire l'optimisation (voir la diapos "Résultats BFGS", cours 4.).
Utilisez la commande hermite.h.quadrature pour calculer $U'(\omega) = E[u'(c)(e^Z-R)]$.
Vous pouvez utiliser le code de mes diapos.

\pagebreak

## Cours 8. La génération de variables aléatoires univariées

### Exercices préliminaires

Si vous avez besoin d'aide pour ces exercices, je suggère le site Wikipédia [ici](https://en.wikipedia.org/wiki/Inverse_transform_sampling).
Sous "References" il y a un lien au livre célèbre et gratuit de [Devroye](http://www.eirene.de/Devroye.pdf) (1986) qui donne plus de détail 

1. Soit $U \sim U([0,1])$ une variable aléatoire uniforme sur l'intervalle $[0,1]$.
Trouvez la densité de la variable aléatoire $X \equiv -(\log U)/\lambda$.
Je suggère les étapes suivantes :
    a. Exprimez $F_X(x)$, la fonction de répartition de $X$ évaluée à $x$, comme la probabilité de l'évènement $X \leq x$.
    a. Replacez $X$ par sa définition, une fonction de $U$.
    a. Exprimez la probabilité en termes de la fonction de répartition de $U$.
    a. Notez que la fonction de répartition d'une v.a. $U([0,1])$ est
    \[
      F_U(u) = \begin{cases} 0 & u \leq 0 \\ u & 0 < u \leq 1 \\ 1 & u > 1 \end{cases}
    \]
    a. Prenez la dérivée de $F_X(x)$ pour obtenir la densité.

1. Supposez que $F$ est une fonction de répartition univariée et que $U \sim U([0,1])$.
Montrez que la fonction de répartition de la variable aléatoire $F^{-1}(U)$ est $F$.

### Exercices théoriques

1. On veut tirer les variables aléatoires de la loi $\mathrm{Ga}(\alpha, 1)$, pour $\alpha \in (0,1)$.
La densité cible est
\[
  f(x) = \frac{1}{\Gamma(\alpha)} x^{\alpha-1} e^{-x}.
\]
On utilise la densité de proposition $\pi g_1(x) + (1-\pi) g_2(x)$, où $\pi \in (0,1)$,
\[
  g_1(x) = \alpha x^{\alpha - 1} 1_{[0,1]}(x), \quad g_2(x) = e^{x-1} 1_{(1,\infty)}(x).
\]
    a. Démontrez que la densité $f(x)$ n'est pas bornée.
    a. Démontrez que la densité $g(x)$ est non-négative et que $\int_0^\infty g(x)\, dx = 1$ pour tous $\pi \in (0,1)$.
    a. Trouvez
    \[
      M_1(\pi) \equiv \sup_{x \in [0,1]} \frac{f(x)}{g(x)}, \quad M_2(\pi) \equiv \sup_{x \in (1,\infty)} \frac{f(x)}{g(x)}.
    \]
    a. Trouvez $\pi \in (0,1)$ qui minimise $M(\pi) = \max(M_1(\pi), M_2(\pi))$.
    a. Quelle est la probabilité d'acceptation si on utilise $g(x)$ (avec $\pi$ optimal) comme densité de proposition pour tirer de la loi $\mathrm{Ga}(\alpha, 1)$, $\alpha \in (0,1)$.
Donnez les chiffres pour $\alpha = 0.01, 0.1, 0.5, 0.9$. La fonction R \texttt{gamma} calcule la fonction $\Gamma(x)$.
    a. Comment peut-on tirer de la loi de proposition (avec densité $g(x)$)?

1. Démontrez que si $X \sim \mathrm{Ga}(\alpha, 1)$, $X/\beta \sim \mathrm{Ga}(\alpha, \beta)$.
Quelle est la pertinence de ce résultat pour le tirage des variables aléatoires gamma.

### Exercice numérique

1. Calculez une approximation au ratio de volumes entre une $d$-hypersphère et un $d$-hypercube.
    a. Exprimez la volume d'une $d$-hypersphère de rayon $\frac{1}{2}$, centrée au point $(\tfrac{1}{2},\ldots,\tfrac{1}{2})$, comme un intégral sur l'hypercube $[0,1]^d$.
    a. Pour $d \in \{4, 8\}$ et $n \in \{1000, 10000\}$, calculez une approximation de l'intégral
        i. avec $n$ $d$-vecteurs pseudo-aléatoires.
        i. avec $n$ $d$-vecteurs quasi-aléatoires du type Halton.
    a. Pour les séquences pseudo-aléatoires, calculez une estimation de l'écart-type de l'approximation.
    a. Comparez vos résultats avec les résultats analytiques [ici](https://en.wikipedia.org/wiki/Volume_of_an_n-ball).

\pagebreak

## Cours 9. La génération de variables aléatoires multivariées

### Exercices préliminaires

Pas d'exercice préliminaire cette semaine, pas de quiz en classe.

### Exercice théorique

1. Pour le modèle de régression gaussien suivant, dérivez les distributions
conditionnelles a posteriori $\beta|h,y$ et $h|\beta,y$.
La loi des données est
\[
  y = X\beta + \epsilon, \quad \epsilon|X \sim N(0,h^{-1}I_n),
\]
et les lois a priori sont
\[
  \beta \sim N(\bar{\beta}, \bar{H}^{-1}_\beta),
  \quad \bar{s}^2h \sim \chi^2(\bar\nu).
\]

### Exercice numérique

1. Avec l'échantillage de Gibbs, simulez la loi a posterior pour
le modèle 3.1 de Wooldridge (2016) "Introductory Econometrics".
(Regardez le fichier \texttt{GPA.R})
Utilisez $M=10000$ tirages, $\bar{\beta} = 0$, $\bar{H}_\beta = 0.01I$, $\bar{s}^2 = 0.01$, $\bar{\nu} = 2$.
Donnez la moyenne et l'écart-type a posteriori des éléments de $\beta$ et les écarts-types numériques pour la moyenne (avec \texttt{mcse} du paquet R \texttt{mcmcse}).

\pagebreak

## Cours 10. La simulation Monte Carlo séquentielle

### Exercice préliminaire

1. Lisez la section "Importance Sampling", pages 293-294 de Judd.
Supposez qu'on veut évaluer l'intégral $E_f[g(x)] = \int f(x) g(x)\, dx$, où $f(x)$
est une densité.
On peut évaluer une fonction $k(x)$ qui est proportionnelle à $f(x)$ : $f(x) = ck(x)$, où $c$ est une constante de normalisation inconnue.
On peut aussi évaluer et simuler une densité normalisée $p(x)$, similaire à la densité $f(x)$, avec le même support et avec la propriété que $f(x)/p(x)$ est borné.
    a. Montrez que $E_p[h_1(x)] = 1/c$ où $h_1(x) = k(x)/p(x)$. Avec un échantillon iid $\{x^{(m)}\}_{m=1}^M$ de la loi avec densité $p(x)$, comment peut-on estimer $1/c$ par Monte Carlo?
    a. Trouvez une fonction $h_2(x)$ telle que $E_p[h_2(x)] = E_f[g(x)]/c$.
    Notez que $E_f[g(x)] = E_p[h_2(x)]/E_p[h_1(x)]$. Avec le même échantillon, comment peut-on estimer $E_f[g(x)]$ comme un ratio d'estimateurs Monte Carlo.

### Autres exercices

Il n'y a pas d'autres exercices pour ce cours.

\pagebreak

## Cours 12, 13. La programmation dynamique

### Exercice numérique

1. $k$ ans après la plantation d'un arbre, la quantité de bois vendable est
\[
  g(k) = e^{-\lambda} \sum_{j=0}^k \frac{\lambda^j}{j!}, \quad \lambda = 30.
\]
On peut utiliser la fonction R \texttt{ppois(k, lambda)} pour évaluer $g(k)$, mais notez qu'il n'y a pas d'interprétation stochastique; la forme sigmoïde de la fonction de répartition d'une variable aléatoire poissonnienne est qualitativement réaliste pour le croissance d'une arbre.
Un propriétaire possède une ferme avec une capacité d'un arbre.
À chaque période $t$, il peut ou couper l'arbre, vendre tout le bois à un prix réel de 1 et planter un nouvel arbre ou laisser continuer à pousser l'arbre. Un arbre planté à $t$ a l'âge 1 à $t+1$.Le taux d'intérêt réel est de $R = 5\%$.
Le propriétaire maximise la valeur actuelle de la récolte du bois.
Il commence avec un arbre d'âge $k_0 = 20$.
    a. Écrivez le problème dynamique du propriétaire dans la forme exposée dans la diapo "La version à horizon infini".
    a. Avec itération par fonction de valeur, trouver la fonction de valeur et la fonction de politique, deux fonctions sur $K \equiv \{1,\ldots,100\}$.