Assignment5.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementing the DGP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100\n",
    "T = 6  # T = 5 in the assignment but since we also included the index 0, we actually have T = 6\n",
    "rho = 0.5\n",
    "pi = 1\n",
    "beta = 1\n",
    "theta = 0\n",
    "np.random.seed(2)\n",
    "def DGP (N, T, rho, pi, beta, theta):\n",
    "    alpha = np.random.normal(size = N)\n",
    "    zeta = np.random.normal(size = N)\n",
    "    nu = np.random.normal(size = N)\n",
    "    eps = np.random.normal(size = N)\n",
    "    x0 = pi*alpha /(1-rho) + (theta*nu + zeta)/np.sqrt(1-rho**2)\n",
    "    y0 = beta*x0 + alpha + eps\n",
    "    x = np.array([x0])\n",
    "    y = np.array([y0])\n",
    "    for i in range(1,T):\n",
    "        xi = np.random.normal(size = N)\n",
    "        xt = rho*x[i-1] + pi*alpha + theta*eps +xi\n",
    "        x = np.append(x, [xt],axis=0)\n",
    "        eps = np.random.normal(size = N)\n",
    "        yt = beta*xt + alpha + eps\n",
    "        y = np.append(y, [yt], axis=0)\n",
    "    x = x.transpose()\n",
    "    y = y.transpose()\n",
    "    return(y, x)\n",
    "\n",
    "def PooledOLS(y, x):\n",
    "    y_pool = y.ravel()\n",
    "    x_pool = x.ravel()\n",
    "    model = sm.OLS(y_pool, x_pool)\n",
    "    results = model.fit()\n",
    "    return([results,x_pool,y_pool])\n",
    "\n",
    "def FD(y, x):\n",
    "    dy = np.diff(y)\n",
    "    dx = np.diff(x)\n",
    "    dy_pool = dy.ravel()\n",
    "    dx_pool = dx.ravel()\n",
    "    model = sm.OLS(dy_pool, dx_pool)\n",
    "    results = model.fit()\n",
    "    return([results, dx_pool, dy_pool])\n",
    "\n",
    "def FE(y, x):\n",
    "    y_FE = (y - y.mean(axis=1, keepdims=True)).ravel()\n",
    "    x_FE = (x - x.mean(axis=1, keepdims=True)).ravel()\n",
    "    model = sm.OLS(y_FE, x_FE)\n",
    "    results = model.fit()\n",
    "    return([results, x_FE, y_FE])\n",
    "\n",
    "def GmmIV(y, x, Z):   # We implemented the IV estimator based on the page 745 of Cameron & Trivedi, because we could not find\n",
    "    y = np.delete(y,0,1)              # the proper package.\n",
    "    x = np.delete(x,0,1)\n",
    "    y = np.diff(y)\n",
    "    x = np.diff(x)\n",
    "    XZ, ZX, ZY, W = 0, 0, 0, 0\n",
    "    T = len(y[0])\n",
    "    N = len(y)\n",
    "    D = np.zeros([T,T+1])\n",
    "    for i in range(T):\n",
    "        D[i][i] = -1\n",
    "        D[i][i+1] = 1\n",
    "    for i in range(N):\n",
    "        XZ = XZ + x[i].dot(Z[i])\n",
    "        ZX = ZX + Z[i].transpose().dot(x[i])\n",
    "        ZY = ZY + Z[i].transpose().dot(y[i])\n",
    "        W = W + Z[i].transpose().dot(D).dot(D.transpose()).dot(Z[i])\n",
    "    if Z[0].shape == (len(Z[0]),):   # some if and else statements are needed because .dot and linalg does not work with scalar \n",
    "        W = 1 / W\n",
    "        Inv = XZ * W * ZX\n",
    "        Inv = 1/Inv\n",
    "        Q = XZ * W * ZY\n",
    "        IV = Inv * Q\n",
    "    else:\n",
    "        W = np.linalg.inv(W)\n",
    "        Inv = XZ.dot(W).dot(ZX)\n",
    "        if Inv.shape == ():\n",
    "            Inv = 1/Inv\n",
    "            Q = XZ.dot(W).dot(ZY)\n",
    "            IV = Inv * Q\n",
    "        else:\n",
    "            Inv = np.linalg.inv(Inv)\n",
    "            Q = XZ.dot(W).dot(ZY)\n",
    "            IV = Inv.dot(Q)\n",
    "    return([IV,x,y])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Exercise 1\n",
    "We first performed Pooled OLS (POLS), FD-OLS and FE-OLS with 1000 simulations. We expect that Pooled OLS is inconsistent, since $x_{it}$ and $\\alpha_i$ are correlated. Therefore we also expect that the POLS to be biased and the t-test is also biased (with a biased t-test of an estimator we mean that the t-test rejects the null-hypothesis H0 : $\\beta = 1$, and hence it is not a good test).\\\n",
    "We do expect that FD-OLS and FE_OLS are consistent and not biased, since these estimators eliminated $\\alpha_i$. Furthermore, we have $\\theta=0$. Therefore, the regressors is not correlated with the error.\\\n",
    "The simulations below confirmed our expectations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 1000\n",
    "beta_est = np.zeros([3,M])\n",
    "beta_cov = np.zeros([3,M])\n",
    "accepted = np.zeros(3)\n",
    "for n in range(M):\n",
    "    [y, x] = DGP(N, T, rho, pi, beta, theta)\n",
    "    result = PooledOLS(y, x)[0]\n",
    "    beta_est[0][n] = result.params[0]\n",
    "    beta_cov[0][n] = result.bse[0]\n",
    "    if result.t_test((1,1)).conf_int()[0][0] <= 1 and result.t_test((1,1)).conf_int()[0][1] >= 1:\n",
    "        accepted[0] += 1\n",
    "    result = FD(y, x)[0]\n",
    "    beta_est[1][n] = result.params[0]\n",
    "    beta_cov[1][n] = result.bse[0]\n",
    "    if result.t_test((1,1)).conf_int()[0][0] <= 1 and result.t_test((1,1)).conf_int()[0][1] >= 1:\n",
    "        accepted[1] += 1\n",
    "    result = FE(y, x)[0]\n",
    "    beta_est[2][n] = result.params[0]\n",
    "    beta_cov[2][n] = result.bse[0]\n",
    "    if result.t_test((1,1)).conf_int()[0][0] <= 1 and result.t_test((1,1)).conf_int()[0][1] >= 1:\n",
    "        accepted[2] += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean of the estimators of respectively POLS, FD and FE:\n",
      "[1.37506661 0.99775135 0.99856202]\n",
      "The average of the standard deviation of the estimators of respectively POLS, FD and FE:\n",
      "[0.01982971 0.05481979 0.04140601]\n",
      "The frequency of not rejecting the nul-hypothesis H0: beta = 1 by the t-test of POLS, FD and FE\n",
      "[0.    0.914 0.92 ]\n"
     ]
    }
   ],
   "source": [
    "print('The mean of the estimators of respectively POLS, FD and FE:')\n",
    "print(np.mean(beta_est, axis = 1))\n",
    "print('The average of the standard deviation of the estimators of respectively POLS, FD and FE:')\n",
    "print(np.mean(beta_cov, axis = 1))\n",
    "print('The frequency of not rejecting the nul-hypothesis H0: beta = 1 by the t-test of POLS, FD and FE')\n",
    "print(accepted/1000)\n",
    "#bins1 = np.histogram(beta_est[0])[1]\n",
    "#plt.hist(beta_est[0], bins = bins1, density = True, histtype = 'bar')\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 2\n",
    "\n",
    "To test the robustness of FD and FE-estimators, we chose some specific values of $\\pi$ that ranges from 0.001 to 1000. We expect that FD and FE are robust to changes of $\\pi$, since they both eliminated $\\alpha_i$ when differencing. The simulation below indeed shows that both estimators are still unbiased.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 1000\n",
    "pi_values = np.array([0.001, 0.1, 5, 50, 1000]) \n",
    "beta_FD = np.zeros([len(pi_values), M])\n",
    "beta_FDsd = np.zeros([len(pi_values), M])\n",
    "beta_FE = np.zeros([len(pi_values), M])\n",
    "beta_FEsd = np.zeros([len(pi_values), M])\n",
    "i = 0\n",
    "for p in pi_values:\n",
    "    for n in range(M):\n",
    "        [y, x] = DGP(N, T, rho, p, beta, theta)\n",
    "        result = FD(y, x)[0]\n",
    "        beta_FD[i][n] = result.params[0]\n",
    "        beta_FDsd[i][n] = result.bse[0]\n",
    "        result = FE(y, x)[0]\n",
    "        beta_FE[i][n] = result.params[0]\n",
    "        beta_FEsd[i][n] = result.bse[0]\n",
    "    i += 1\n",
    "       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean of the estimator of FD over 5 values of pi:\n",
      "[[1.00155294]\n",
      " [1.00190575]\n",
      " [0.99938626]\n",
      " [1.00182775]\n",
      " [0.99989678]]\n",
      "The average of the standard deviation of the estimator FD over 5 values of pi:\n",
      "[[0.05482207]\n",
      " [0.05488792]\n",
      " [0.05479792]\n",
      " [0.05485425]\n",
      " [0.05490479]]\n",
      "The mean of the estimator of FE over 5 values of pi:\n",
      "[[1.0010828 ]\n",
      " [1.00087461]\n",
      " [0.99850841]\n",
      " [1.00132603]\n",
      " [1.00103798]]\n",
      "The average of the standard deviation of the estimator FE over 5 values of pi:\n",
      "[[0.04135248]\n",
      " [0.04144228]\n",
      " [0.04142116]\n",
      " [0.04153766]\n",
      " [0.04147161]]\n"
     ]
    }
   ],
   "source": [
    "print('The mean of the estimator of FD over 5 values of pi:')\n",
    "print(np.mean(beta_FD, axis = 1, keepdims = True))\n",
    "print('The average of the standard deviation of the estimator FD over 5 values of pi:')\n",
    "print(np.mean(beta_FDsd, axis = 1, keepdims = True))\n",
    "print('The mean of the estimator of FE over 5 values of pi:')\n",
    "print(np.mean(beta_FE, axis = 1, keepdims = True))\n",
    "print('The average of the standard deviation of the estimator FE over 5 values of pi:')\n",
    "print(np.mean(beta_FEsd, axis = 1, keepdims = True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 3\n",
    "We analyzed the impact of $\\theta = 0.5$ for the FD and FE estimators. Since $\\theta\\neq 0$, we now have that $x_{it}$ is correlated with $\\epsilon_{it-1}$. Therefore $\\Delta x_{it}$ is correlated with $\\Delta \\epsilon_{it}$ and $x_i - 1/T \\sum_t x_{it}$ is correlated with $\\epsilon_i -1/T\\sum_t \\epsilon_{it}$. Therefore, we expect that FD and FE estimators are biased. The results of the simulation below agree with our expectation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 1000 \n",
    "beta_FD = np.zeros(M)\n",
    "beta_FDsd = np.zeros(M)\n",
    "beta_FE = np.zeros(M)\n",
    "beta_FEsd = np.zeros(M)\n",
    "theta = 0.5\n",
    "for n in range(M):\n",
    "    [y, x] = DGP(N, T, rho, pi, beta, theta)\n",
    "    result = FD(y, x)[0]\n",
    "    beta_FD[n] = result.params[0]\n",
    "    beta_FDsd[n] = result.bse[0]\n",
    "    result = FE(y, x)[0]\n",
    "    beta_FE[n] = result.params[0]\n",
    "    beta_FEsd[n] = result.bse[0]    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean of the estimator of FD for theta = 0.5:\n",
      "0.7012540872881287\n",
      "The average of the standard deviation of the estimator FD:\n",
      "0.047127249694586405\n",
      "The mean of the estimator of FE for theta = 5:\n",
      "0.8896116765867669\n",
      "The average of the standard deviation of the estimator FE:\n",
      "0.036660910060674894\n"
     ]
    }
   ],
   "source": [
    "print('The mean of the estimator of FD for theta = 0.5:')\n",
    "print(np.mean(beta_FD))\n",
    "print('The average of the standard deviation of the estimator FD:')\n",
    "print(np.mean(beta_FDsd))\n",
    "print('The mean of the estimator of FE for theta = 5:')\n",
    "print(np.mean(beta_FE))\n",
    "print('The average of the standard deviation of the estimator FE:')\n",
    "print(np.mean(beta_FEsd))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 4\n",
    "We impletemented and printed the $Z_i$ matrices. We are considering three instrument matrices which we denote $Z_1$, $Z_2$ and $Z_3$. They are matrices that respectively uses 1, $(T-1)$ and $1/2T(T-1)$ instruments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Z-matrix when we use 1 instrument:\n",
      "[[ 2.21535107 -0.01557787 -0.77228213 -0.23094584]\n",
      " [ 1.63046986  3.60999553  0.79243676  3.80038479]\n",
      " [ 2.00026584  2.96840115  2.89115519  3.55719188]\n",
      " [-1.56789146 -1.28296065  0.25210545  1.02508187]\n",
      " [ 1.21798656  1.91053711  2.22455612  2.37959091]\n",
      " [-0.34650266  1.08739286  0.65036611  2.8645199 ]\n",
      " [ 2.40417779  3.14835973  3.53576829  2.85063266]\n",
      " [ 1.35072159 -0.47551072  2.53208371  2.47842111]\n",
      " [-1.05250117 -0.9149518   1.58669109 -0.68628944]\n",
      " [-3.53983405 -5.14037862 -3.37118937 -5.73196952]\n",
      " [ 2.86860506  4.79682918  3.29959592  4.90200808]\n",
      " [ 0.31612718 -0.67670241  0.16830319 -1.46905925]\n",
      " [ 2.91277602  3.64572824  5.37423391  3.86224645]\n",
      " [ 2.28590472  3.01253398  1.11498655  0.3309319 ]\n",
      " [-1.80797931 -4.73875971 -1.14360492 -0.01423225]\n",
      " [ 2.29316251  0.46358194  1.26260989  1.72560399]\n",
      " [-2.69320628 -1.55791531 -0.50793143 -1.61538   ]\n",
      " [-1.79028859 -3.14514583  0.06798962  0.29246337]\n",
      " [ 2.66284266  2.70419502  1.62678984  3.25530067]\n",
      " [ 3.06175375  1.5540789   1.33368964 -0.56258721]\n",
      " [ 2.45731309  2.57880621  2.81799226  0.67552297]\n",
      " [ 1.77988382  1.90363519  0.88944051  1.75567182]\n",
      " [ 3.18405813  4.13496122  2.56884591  2.74222504]\n",
      " [ 0.4362901   1.26297319  2.43631205  3.21755629]\n",
      " [-2.75923016 -0.31746152  0.3680428  -0.29380835]\n",
      " [-3.77877499 -2.88902554 -5.57535159 -6.79344279]\n",
      " [-1.20066208  0.03034406  0.88536295  0.71196933]\n",
      " [-1.48005706 -1.62838743 -1.38677856 -0.45558981]\n",
      " [-3.35174006 -2.11081855 -4.15440451 -2.38032623]\n",
      " [-3.04825266 -3.27714658 -4.6964954  -3.38570032]\n",
      " [ 3.77289231  2.36754006  2.98991475  3.8783681 ]\n",
      " [-3.69107727 -3.77714854 -4.68804725 -5.27080053]\n",
      " [ 0.23875501 -1.64903617 -2.58973078 -4.98409807]\n",
      " [ 0.70892736 -0.22090159  1.35602281  0.93085839]\n",
      " [-0.14856185 -1.60177682 -0.43949107  0.16013518]\n",
      " [-1.32084288 -3.1111166  -2.384272   -1.56222454]\n",
      " [-2.70098925 -2.80245614 -3.98861947 -4.37118651]\n",
      " [-2.44740914 -4.4619203  -3.71324071 -2.17525332]\n",
      " [ 0.27412525 -0.06295749 -1.71346081  0.32495598]\n",
      " [-2.37258559 -0.76286959 -0.44613863  0.08045798]\n",
      " [-0.57355524 -0.84745338 -0.41416013 -0.11442039]\n",
      " [ 2.4792351   2.24624435  1.99819614  2.09541119]\n",
      " [ 4.38946657  5.40688719  4.53758097  3.71846438]\n",
      " [ 5.23505511  4.46840062  2.73491678  3.88654339]\n",
      " [ 1.13368671  1.8439435  -0.58042388 -1.51135227]\n",
      " [ 2.0308354   2.86053724  2.13595993  1.77591766]\n",
      " [-3.47797751 -5.46480544 -4.62652958 -6.76723109]\n",
      " [-3.96844449 -4.79285518 -5.10237774 -3.84043717]\n",
      " [-1.31671147  0.37990497  1.75422736  0.44353381]\n",
      " [ 0.72252336  0.49334156  1.36494946  1.82732964]\n",
      " [ 2.85487064  2.73132508 -0.1130353  -0.25806287]\n",
      " [ 2.09854205  3.98729388  2.74627042  1.82999291]\n",
      " [-2.2500852  -0.71916972 -0.34705676 -1.22205858]\n",
      " [-1.33272967 -2.49677786 -1.56435752 -0.84525684]\n",
      " [ 6.36383987  5.92127584  5.77995504  5.31028202]\n",
      " [-1.37703686 -1.39484348 -2.53236648 -0.84889468]\n",
      " [-1.20126433 -2.20318338 -2.58533935 -4.610308  ]\n",
      " [ 1.36537247  3.03661302  1.98157414  1.57745155]\n",
      " [ 1.35352868 -0.14823647  1.00045163  2.10943826]\n",
      " [-1.79677936 -3.18292377 -0.45225152 -1.32722784]\n",
      " [ 3.1306292   0.94886208  2.58464082  3.05521344]\n",
      " [-4.75114693 -4.33212399 -3.80849538 -2.74775199]\n",
      " [ 0.43142439  2.27835197  1.13135005  0.34421854]\n",
      " [-0.75564153 -1.97572884 -0.96733322 -2.00333994]\n",
      " [-0.08826621 -0.53360072 -1.04783515 -2.9284048 ]\n",
      " [ 2.96172865  2.44620677  3.13630463  3.15990156]\n",
      " [-0.94008689 -1.17119725 -1.88232577 -1.89243466]\n",
      " [-3.74135646 -2.52845393 -1.32643566  0.31106371]\n",
      " [ 1.67213284  2.54057029  1.44333826  0.71256039]\n",
      " [ 0.17408577 -2.32551161 -1.03786768 -1.51183256]\n",
      " [-2.13847935 -4.79012408 -4.73920109 -1.75198555]\n",
      " [-1.9598903  -0.16086293  0.25188577 -1.18563605]\n",
      " [-1.36132142  1.20270554  0.96988824  2.15155726]\n",
      " [-4.18494703 -3.25910955 -4.59417388 -2.65070213]\n",
      " [ 0.23868156 -1.2879731   0.12482431  0.81769239]\n",
      " [ 1.40411112  3.15014573  0.4886103   1.79836016]\n",
      " [-2.06875792 -2.50709348 -0.85005405 -1.96787205]\n",
      " [-3.39779609 -0.95643279 -2.07017758 -3.76160647]\n",
      " [-2.24169174 -1.02648461  0.92038755  1.25314059]\n",
      " [-0.17972131  0.50849619 -1.39096275 -0.9620043 ]\n",
      " [-2.1120118   0.04763135 -0.38222747 -0.04375521]\n",
      " [ 1.97591232  2.25374456  2.45419831  2.14911965]\n",
      " [ 0.41720453  1.44061358  0.2078086   0.48630838]\n",
      " [ 3.61056029  2.23566086  2.11902602 -0.98132822]\n",
      " [ 0.54513639  0.08458025 -0.20306213  1.52832078]\n",
      " [ 4.3395511   2.74639511  1.94877495  2.87824104]\n",
      " [-1.76605542 -1.15254423 -2.33810644 -2.82802582]\n",
      " [ 2.4005693   0.77216017  3.65129735  2.87590004]\n",
      " [-2.08042808 -2.35137581 -0.40625958 -1.16161929]\n",
      " [-1.38001304 -0.98580674 -0.24716385  0.4376449 ]\n",
      " [ 0.37176041  0.24111125 -1.20455619 -0.44442009]\n",
      " [ 3.03092335  3.8250164   4.82109888  2.21902027]\n",
      " [-1.35095902 -1.74789795 -1.94171453 -4.38629671]\n",
      " [ 4.85415802  4.04771855  3.29911606  3.404481  ]\n",
      " [-0.73466256  0.40286174  0.4585051  -0.67043719]\n",
      " [-3.32428809 -2.59869433 -0.44194006 -4.76798996]\n",
      " [ 2.11737688  3.57362907  0.88608557  1.72098353]\n",
      " [ 0.98186284  1.1084157  -0.06411317  0.51173003]\n",
      " [-0.39880764 -0.85305605 -0.25585658 -0.80254863]\n",
      " [-2.1670597  -0.78986448 -1.36421745 -1.3701811 ]]\n",
      "The Z-matrix when we use T-1 instruments:\n",
      "[[[ 2.21535107  0.          0.          0.        ]\n",
      "  [ 0.         -0.01557787  0.          0.        ]\n",
      "  [ 0.          0.         -0.77228213  0.        ]\n",
      "  [ 0.          0.          0.         -0.23094584]]\n",
      "\n",
      " [[ 1.63046986  0.          0.          0.        ]\n",
      "  [ 0.          3.60999553  0.          0.        ]\n",
      "  [ 0.          0.          0.79243676  0.        ]\n",
      "  [ 0.          0.          0.          3.80038479]]\n",
      "\n",
      " [[ 2.00026584  0.          0.          0.        ]\n",
      "  [ 0.          2.96840115  0.          0.        ]\n",
      "  [ 0.          0.          2.89115519  0.        ]\n",
      "  [ 0.          0.          0.          3.55719188]]\n",
      "\n",
      " ...\n",
      "\n",
      " [[ 0.98186284  0.          0.          0.        ]\n",
      "  [ 0.          1.1084157   0.          0.        ]\n",
      "  [ 0.          0.         -0.06411317  0.        ]\n",
      "  [ 0.          0.          0.          0.51173003]]\n",
      "\n",
      " [[-0.39880764  0.          0.          0.        ]\n",
      "  [ 0.         -0.85305605  0.          0.        ]\n",
      "  [ 0.          0.         -0.25585658  0.        ]\n",
      "  [ 0.          0.          0.         -0.80254863]]\n",
      "\n",
      " [[-2.1670597   0.          0.          0.        ]\n",
      "  [ 0.         -0.78986448  0.          0.        ]\n",
      "  [ 0.          0.         -1.36421745  0.        ]\n",
      "  [ 0.          0.          0.         -1.3701811 ]]]\n",
      "The Z-matrix when we use 1/2T(T-1) instruments:\n",
      "[[[ 2.21535107  0.          0.         ...  0.          0.\n",
      "    0.        ]\n",
      "  [ 0.          2.21535107 -0.01557787 ...  0.          0.\n",
      "    0.        ]\n",
      "  [ 0.          0.          2.21535107 ... -0.77228213  0.\n",
      "    0.        ]\n",
      "  [ 0.          0.          0.         ... -0.01557787 -0.77228213\n",
      "   -0.23094584]]\n",
      "\n",
      " [[ 1.63046986  0.          0.         ...  0.          0.\n",
      "    0.        ]\n",
      "  [ 0.          1.63046986  3.60999553 ...  0.          0.\n",
      "    0.        ]\n",
      "  [ 0.          0.          1.63046986 ...  0.79243676  0.\n",
      "    0.        ]\n",
      "  [ 0.          0.          0.         ...  3.60999553  0.79243676\n",
      "    3.80038479]]\n",
      "\n",
      " [[ 2.00026584  0.          0.         ...  0.          0.\n",
      "    0.        ]\n",
      "  [ 0.          2.00026584  2.96840115 ...  0.          0.\n",
      "    0.        ]\n",
      "  [ 0.          0.          2.00026584 ...  2.89115519  0.\n",
      "    0.        ]\n",
      "  [ 0.          0.          0.         ...  2.96840115  2.89115519\n",
      "    3.55719188]]\n",
      "\n",
      " ...\n",
      "\n",
      " [[ 0.98186284  0.          0.         ...  0.          0.\n",
      "    0.        ]\n",
      "  [ 0.          0.98186284  1.1084157  ...  0.          0.\n",
      "    0.        ]\n",
      "  [ 0.          0.          0.98186284 ... -0.06411317  0.\n",
      "    0.        ]\n",
      "  [ 0.          0.          0.         ...  1.1084157  -0.06411317\n",
      "    0.51173003]]\n",
      "\n",
      " [[-0.39880764  0.          0.         ...  0.          0.\n",
      "    0.        ]\n",
      "  [ 0.         -0.39880764 -0.85305605 ...  0.          0.\n",
      "    0.        ]\n",
      "  [ 0.          0.         -0.39880764 ... -0.25585658  0.\n",
      "    0.        ]\n",
      "  [ 0.          0.          0.         ... -0.85305605 -0.25585658\n",
      "   -0.80254863]]\n",
      "\n",
      " [[-2.1670597   0.          0.         ...  0.          0.\n",
      "    0.        ]\n",
      "  [ 0.         -2.1670597  -0.78986448 ...  0.          0.\n",
      "    0.        ]\n",
      "  [ 0.          0.         -2.1670597  ... -1.36421745  0.\n",
      "    0.        ]\n",
      "  [ 0.          0.          0.         ... -0.78986448 -1.36421745\n",
      "   -1.3701811 ]]]\n"
     ]
    }
   ],
   "source": [
    "theta = 0.5\n",
    "[y, x] = DGP(N, T, rho, pi, beta, theta)\n",
    "Z1 = np.zeros([N,T-2])\n",
    "for i in range(N):\n",
    "    for t in range((T-2)):\n",
    "        Z1[i][t] = x[i][t+1]\n",
    "print('The Z-matrix when we use 1 instrument:')\n",
    "print(Z1)\n",
    "Z2 = np.zeros([N, T-2, T-2])\n",
    "for i in range(N):\n",
    "    for t in range((T-2)):\n",
    "        Z2[i][t][t] = x[i][t+1]\n",
    "r = 1 + (T-3) * 2\n",
    "print('The Z-matrix when we use T-1 instruments:')\n",
    "print(Z2)\n",
    "Z3 = np.zeros([N, T-2, r])\n",
    "for i in range(N):\n",
    "    Z3[i][0][0] = x[i][1]\n",
    "    for t in range(1,(T-2)):\n",
    "        Z3[i][t][t:(2*t+1)] = x[i][1:(t+2)]\n",
    "print('The Z-matrix when we use 1/2T(T-1) instruments:')\n",
    "print(Z3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Exercise 5\n",
    "We looked at the bias and standard-deviation of the IV-estimator using instrument matrix Z1 (with 1 instrument), Z2 (with T-1 instruments) and Z3 (with $\\frac{1}{2}T(T-1)$ instruments). We expect that these are all unbiased, since our instrument is uncorrelated with the errors. We have also varied $\\theta$ and we observed that all three estimators are unbiased for a range of $\\theta$ from 0.001 to 1000. We did not observe any preference for the estimators, in terms of unbiasedness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 1000\n",
    "thetas = [0.001, 0.5, 20, 1000]\n",
    "beta_IV = np.zeros([len(thetas), 3, M])\n",
    "j = 0\n",
    "r = 1 + (T-3) * 2\n",
    "Z1 = np.zeros([N,T-2])\n",
    "Z2 = np.zeros([N, T-2, T-2])\n",
    "Z3 = np.zeros([N, T-2, r])\n",
    "for theta in thetas:\n",
    "    for n in range(M):\n",
    "        [y, x] = DGP(N, T, rho, pi, beta, theta)\n",
    "        for i in range(N):\n",
    "            for t in range((T-2)):\n",
    "                Z1[i][t] = x[i][t+1]\n",
    "                Z2[i][t][t] = x[i][t+1]\n",
    "        for i in range(N):\n",
    "            Z3[i][0][0] = x[i][1]\n",
    "            for t in range(1,(T-2)):\n",
    "                Z3[i][t][t:(2*t+1)] = x[i][1:(t+2)]\n",
    "        beta_IV[j][0][n] = GmmIV(y, x, Z1)[0]\n",
    "        beta_IV[j][1][n] = GmmIV(y, x, Z2)[0]\n",
    "        beta_IV[j][2][n] = GmmIV(y, x, Z3)[0]\n",
    "    j += 1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 0.001\n",
      "[0.99927933 0.99978592 0.9990809 ]\n",
      "The standard deviation of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 0.001\n",
      "[0.14786741 0.13708783 0.15489238]\n",
      "The mean of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 0.5\n",
      "[0.99679396 0.98833223 0.97478881]\n",
      "The standard deviation of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 0.5\n",
      "[0.12189048 0.11478969 0.13046953]\n",
      "The mean of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 20\n",
      "[1.00011015 0.99957288 0.99918854]\n",
      "The standard deviation of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 20\n",
      "[0.00483556 0.00476227 0.00408477]\n",
      "The mean of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 1000\n",
      "[1.00000508 0.99999434 0.99998195]\n",
      "The standard deviation of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 1000\n",
      "[1.00410262e-04 9.80324832e-05 8.11049007e-05]\n"
     ]
    }
   ],
   "source": [
    "print('The mean of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 0.001')\n",
    "print(np.mean(beta_IV[0], axis = 1))\n",
    "print('The standard deviation of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 0.001')\n",
    "print(np.std(beta_IV[0], axis = 1))\n",
    "print('The mean of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 0.5')\n",
    "print(np.mean(beta_IV[1], axis = 1))\n",
    "print('The standard deviation of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 0.5')\n",
    "print(np.std(beta_IV[1], axis = 1))\n",
    "print('The mean of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 20')\n",
    "print(np.mean(beta_IV[2], axis = 1))\n",
    "print('The standard deviation of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 20')\n",
    "print(np.std(beta_IV[2], axis = 1))\n",
    "print('The mean of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 1000')\n",
    "print(np.mean(beta_IV[3], axis = 1))\n",
    "print('The standard deviation of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 1000')\n",
    "print(np.std(beta_IV[3], axis = 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 6\n",
    "We changed the sample size from 100 to 1000 and we performed the above experiment with $\\theta = 0.5$ over 1000 simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 1000\n",
    "theta = 0.5\n",
    "N = 1000\n",
    "beta_IV = np.zeros([3, M])\n",
    "r = 1 + (T-3) * 2\n",
    "Z1 = np.zeros([N,T-2])\n",
    "Z2 = np.zeros([N, T-2, T-2])\n",
    "Z3 = np.zeros([N, T-2, r])\n",
    "\n",
    "for n in range(M):\n",
    "    [y, x] = DGP(N, T, rho, pi, beta, theta)\n",
    "    for i in range(N):\n",
    "        for t in range((T-2)):\n",
    "            Z1[i][t] = x[i][t+1]\n",
    "            Z2[i][t][t] = x[i][t+1]\n",
    "    for i in range(N):\n",
    "        Z3[i][0][0] = x[i][1]\n",
    "        for t in range(1,(T-2)):\n",
    "            Z3[i][t][t:(2*t+1)] = x[i][1:(t+2)]\n",
    "    beta_IV[0][n] = GmmIV(y, x, Z1)[0]\n",
    "    beta_IV[1][n] = GmmIV(y, x, Z2)[0]\n",
    "    beta_IV[2][n] = GmmIV(y, x, Z3)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean (over 1000 samples) of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 0.5\n",
      "[1.00067827 0.99996193 0.9985324 ]\n",
      "The standard deviation (over 1000 samples) of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 0.5\n",
      "[0.03896822 0.03715333 0.0417763 ]\n"
     ]
    }
   ],
   "source": [
    "print('The mean (over 1000 samples) of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 0.5')\n",
    "print(np.mean(beta_IV, axis = 1))\n",
    "print('The standard deviation (over 1000 samples) of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 0.5')\n",
    "print(np.std(beta_IV, axis = 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, the standard deviation of the estimators are getting smaller due to the larger sample size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 7\n",
    "To investigate persistence we set $\\rho = 0.99$. Then $x_{it-1}$ because relatively large with respect to $\\epsilon_{it-1}$, $\\alpha_i$ and $\\xi_{it}$. Since also $\\rho\\to 1$, this makes $x_{it}$ and $x_{it-1}$ similar and thereby induces persistance. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 100\n",
    "theta = 0.5\n",
    "N = 1000\n",
    "rho = 0.99\n",
    "beta_IV = np.zeros([3, M])\n",
    "r = 1 + (T-3) * 2\n",
    "Z1 = np.zeros([N,T-2])\n",
    "Z2 = np.zeros([N, T-2, T-2])\n",
    "Z3 = np.zeros([N, T-2, r])\n",
    "for n in range(M):\n",
    "    [y, x] = DGP(N, T, rho, pi, beta, theta)\n",
    "    for i in range(N):\n",
    "        for t in range((T-2)):\n",
    "            Z1[i][t] = x[i][t+1]\n",
    "            Z2[i][t][t] = x[i][t+1]\n",
    "    for i in range(N):\n",
    "        Z3[i][0][0] = x[i][1]\n",
    "        for t in range(1,(T-2)):\n",
    "            Z3[i][t][t:(2*t+1)] = x[i][1:(t+2)]\n",
    "    beta_IV[0][n] = GmmIV(y, x, Z1)[0]\n",
    "    beta_IV[1][n] = GmmIV(y, x, Z2)[0]\n",
    "    beta_IV[2][n] = GmmIV(y, x, Z3)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean (over 1000 samples) of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 0.5\n",
      "[0.12649626 0.83238908 0.74822692]\n",
      "The standard deviation (over 1000 samples) of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 0.5\n",
      "[4.29736008 0.83362386 0.44803678]\n"
     ]
    }
   ],
   "source": [
    "print('The mean (over 1000 samples) of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 0.5')\n",
    "print(np.mean(beta_IV, axis = 1))\n",
    "print('The standard deviation (over 1000 samples) of the IV-estimators, respectively Z1, Z2 and Z3, calculated for theta = 0.5')\n",
    "print(np.std(beta_IV, axis = 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimators are highly biased. This shows that our instrument is weak when we have persistance. This is because persistance makes $\\Delta x_{it}$ small, which makes estimating $\\beta$ hard."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 8\n",
    "In this exercise we aim to estimate $\\rho$. We regressed on $x_{it-1}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean and the average standard deviation of the POLS estimator over 100 simulations:\n",
      "[0.8703968678639137, 0.021861854771394632]\n",
      "The mean and the average standard deviation of the FD estimator over 100 simulations:\n",
      "[-0.24919912909517816, 0.04836349169687822]\n",
      "The mean and the average standard deviation of the FE estimator over 100 simulations:\n",
      "[0.16981040037700112, 0.04414082378295921]\n"
     ]
    }
   ],
   "source": [
    "M = 100\n",
    "N = 100\n",
    "T = 6  \n",
    "rho = 0.5\n",
    "pi = 1\n",
    "beta = 1\n",
    "theta = 0\n",
    "beta_pool = np.zeros([2,M])\n",
    "beta_FD = np.zeros([2,M])\n",
    "beta_FE = np.zeros([2,M])\n",
    "\n",
    "for n in range(M):\n",
    "    x = DGP(N, T, rho, pi, beta, theta)[1]\n",
    "    Y = np.delete(x, 0,1)\n",
    "    X = np.delete(x, (len(x[0])-1),1)\n",
    "    resultpool = PooledOLS (Y, X)[0]\n",
    "    resultFD = FD(Y, X)[0]\n",
    "    resultFE = FE(Y, X)[0]\n",
    "    beta_pool[0][n] = resultpool.params\n",
    "    beta_pool[1][n] = resultpool.bse\n",
    "    beta_FD[0][n] = resultFD.params\n",
    "    beta_FD[1][n] = resultFD.bse\n",
    "    beta_FE[0][n] = resultFE.params\n",
    "    beta_FE[1][n] = resultFE.bse\n",
    "print(\"The mean and the average standard deviation of the POLS estimator over 100 simulations:\")\n",
    "print([np.mean(beta_pool[0]), np.mean(beta_pool[1])])\n",
    "print(\"The mean and the average standard deviation of the FD estimator over 100 simulations:\")\n",
    "print([np.mean(beta_FD[0]), np.mean(beta_FD[1])])\n",
    "print(\"The mean and the average standard deviation of the FE estimator over 100 simulations:\")\n",
    "print([np.mean(beta_FE[0]), np.mean(beta_FE[1])])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the estimators are highly biased. This is also explained in Cameron and Trivedi. The Pooled OLS is biased because $x_{it-1}$ is correlated with $\\alpha_i$. The FD-estimator is biased because $\\Delta x_{it-1}$ is correlated with $\\Delta \\xi_{it}$ since $x_{it-1}$ is correlated with $\\xi_{it-1}$, which is also the same reason that the FE (within) estimator is biased."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we use a large $T$ then we get a long panel data. It is stated in Cameron and Trivedi that the biasedness of the FE (within) estimator will improve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean and the average standard deviation of the POLS estimator over 100 simulations:\n",
      "[0.8743835171738656, 0.00153267333893189]\n",
      "The mean and the average standard deviation of the FD estimator over 100 simulations:\n",
      "[-0.25015045528972235, 0.003064817154493543]\n",
      "The mean and the average standard deviation of the FE estimator over 100 simulations:\n",
      "[0.4981867826194995, 0.0027432312366061605]\n"
     ]
    }
   ],
   "source": [
    "M = 100\n",
    "N = 100\n",
    "T = 1000 \n",
    "rho = 0.5\n",
    "pi = 1\n",
    "beta = 1\n",
    "theta = 0\n",
    "beta_pool = np.zeros([2,M])\n",
    "beta_FD = np.zeros([2,M])\n",
    "beta_FE = np.zeros([2,M])\n",
    "\n",
    "for n in range(M):\n",
    "    x = DGP(N, T, rho, pi, beta, theta)[1]\n",
    "    Y = np.delete(x, 0,1)\n",
    "    X = np.delete(x, (len(x[0])-1),1)\n",
    "    resultpool = PooledOLS (Y, X)[0]\n",
    "    resultFD = FD(Y, X)[0]\n",
    "    resultFE = FE(Y, X)[0]\n",
    "    beta_pool[0][n] = resultpool.params\n",
    "    beta_pool[1][n] = resultpool.bse\n",
    "    beta_FD[0][n] = resultFD.params\n",
    "    beta_FD[1][n] = resultFD.bse\n",
    "    beta_FE[0][n] = resultFE.params\n",
    "    beta_FE[1][n] = resultFE.bse\n",
    "print(\"The mean and the average standard deviation of the POLS estimator over 100 simulations:\")\n",
    "print([np.mean(beta_pool[0]), np.mean(beta_pool[1])])\n",
    "print(\"The mean and the average standard deviation of the FD estimator over 100 simulations:\")\n",
    "print([np.mean(beta_FD[0]), np.mean(beta_FD[1])])\n",
    "print(\"The mean and the average standard deviation of the FE estimator over 100 simulations:\")\n",
    "print([np.mean(beta_FE[0]), np.mean(beta_FE[1])])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed, we observe that the within estimator has improved, since its estimated $\\rho$ is around $0.498$, which is much closer to the true value $0.5$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we also increased the sample size $N$ to 1000. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean and the average standard deviation of the POLS estimator over 100 simulations:\n",
      "[0.8758107861216111, 0.00682078749479391]\n",
      "The mean and the average standard deviation of the FD estimator over 100 simulations:\n",
      "[-0.2484408030405424, 0.015297953487593594]\n",
      "The mean and the average standard deviation of the FE estimator over 100 simulations:\n",
      "[0.16923850520170988, 0.013922870995855229]\n"
     ]
    }
   ],
   "source": [
    "M = 100\n",
    "N = 1000\n",
    "T = 6 \n",
    "rho = 0.5\n",
    "pi = 1\n",
    "beta = 1\n",
    "theta = 0\n",
    "beta_pool = np.zeros([2,M])\n",
    "beta_FD = np.zeros([2,M])\n",
    "beta_FE = np.zeros([2,M])\n",
    "\n",
    "for n in range(M):\n",
    "    x = DGP(N, T, rho, pi, beta, theta)[1]\n",
    "    Y = np.delete(x, 0,1)\n",
    "    X = np.delete(x, (len(x[0])-1),1)\n",
    "    resultpool = PooledOLS (Y, X)[0]\n",
    "    resultFD = FD(Y, X)[0]\n",
    "    resultFE = FE(Y, X)[0]\n",
    "    beta_pool[0][n] = resultpool.params\n",
    "    beta_pool[1][n] = resultpool.bse\n",
    "    beta_FD[0][n] = resultFD.params\n",
    "    beta_FD[1][n] = resultFD.bse\n",
    "    beta_FE[0][n] = resultFE.params\n",
    "    beta_FE[1][n] = resultFE.bse\n",
    "print(\"The mean and the average standard deviation of the POLS estimator over 100 simulations:\")\n",
    "print([np.mean(beta_pool[0]), np.mean(beta_pool[1])])\n",
    "print(\"The mean and the average standard deviation of the FD estimator over 100 simulations:\")\n",
    "print([np.mean(beta_FD[0]), np.mean(beta_FD[1])])\n",
    "print(\"The mean and the average standard deviation of the FE estimator over 100 simulations:\")\n",
    "print([np.mean(beta_FE[0]), np.mean(beta_FE[1])])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In contrast to increasing $T$, we did not observe any changes in the biasedness."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}