diff --git a/data/US_Data.xlsx b/data/US_Data.xlsx new file mode 100644 index 0000000..0ccbbcf Binary files /dev/null and b/data/US_Data.xlsx differ diff --git a/exercise05-jl.ipynb b/exercise05-jl.ipynb new file mode 100644 index 0000000..944e2cc --- /dev/null +++ b/exercise05-jl.ipynb @@ -0,0 +1,346 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Exercise 5\n", + "- Required packages:\n", + " - `Plots.jl` (for plotting)\n", + " - `ExcelReaders.jl` (to read excel file)" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "using Plots\n", + "pyplot()\n", + "using ExcelReaders" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Functions for Hodrick-Prescott (HP) filter\n", + "This exercise allows us to showcase how different methods of computing the HP filter affects runtime and memory use.\n", + "\n", + "- First, we'll use the loop method with a dense matrix to create matrix $H$.\n", + "- Second, we'll do a loop again, but with a sparse matrix.\n", + "- Third, we'll vectorize the calculation with a dense matrix and via the function `spdiagm` to create matrix $H$.\n", + "- Finally, we'll vectorize the calculation with a sparse matrix and via the function `spdiagm` to create matrix $H$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Declare the matrices" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "abstract type MatrixType end\n", + "struct Dense <: MatrixType end\n", + "struct Sparse <: MatrixType end\n", + "abstract type AssignMethod end\n", + "struct Loop <: AssignMethod end\n", + "struct Vectorized <: AssignMethod end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Loop version" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "function create_matrix(λ, N::Integer, ::Dense, ::Loop)\n", + " H = zeros(N, N)\n", + " assign!(H, λ, N)\n", + " return H\n", + "end;\n", + "\n", + "function create_matrix(λ, N::Integer, ::Sparse, ::Loop)\n", + " H = spzeros(N, N)\n", + " assign!(H, λ, N)\n", + " return H\n", + "end;\n", + "\n", + "function assign!(H::AbstractMatrix, λ, N) \n", + " for j = 1:N\n", + " for i = 1:N\n", + " if j == 1\n", + " H[1, 1] = 1+λ;\n", + " H[2, 1] = -2λ;\n", + " H[3, 1] = λ;\n", + " elseif j == 2\n", + " H[1, 2] = -2λ;\n", + " H[2, 2] = 1 + 5λ;\n", + " H[3, 2] = -4λ;\n", + " H[4, 2] = λ;\n", + " elseif j == N - 1\n", + " H[N - 3, N - 1] = λ;\n", + " H[N - 2, N - 1] = -4λ;\n", + " H[N - 1, N - 1] = 1 + 5λ;\n", + " H[N, N - 1] = -2λ;\n", + " elseif j == N\n", + " H[N - 2, N] = λ;\n", + " H[N - 1, N] = -2λ;\n", + " H[N, N] = 1+λ;\n", + " else\n", + " H[j - 2, j] = λ;\n", + " H[j - 1, j] = -4λ;\n", + " H[j, j] = 1 + 6λ;\n", + " H[j + 1, j] = -4λ;\n", + " H[j + 2, j] = λ;\n", + " end\n", + " end\n", + " end\n", + "end;" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Vectorized version" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "function create_matrix(λ::Real, N::Integer, ::Sparse, ::Vectorized) \n", + " return spdiagm(-2 => fill(λ, N-2),\n", + " -1 => vcat(-2λ, fill(-4λ, N - 3), -2λ),\n", + " 0 => vcat(1 + λ, 1 + 5λ, fill(1 + 6λ, N-4),\n", + " 1 + 5λ, 1 + λ),\n", + " 1 => vcat(-2λ, fill(-4λ, N - 3), -2λ),\n", + " 2 => fill(λ, N-2))\n", + "end;\n", + "\n", + "function create_matrix(λ::Real, N::Integer, ::Dense, ::Vectorized)\n", + " H = zeros(N, N)\n", + " H += diagm(fill(λ, N-2), -2)\n", + " H += diagm(vcat(-2λ, fill(-4λ, N - 3), -2λ), -1)\n", + " H += diagm(vcat(1 + λ, 1 + 5λ, fill(1 + 6λ, N-4),\n", + " 1 + 5λ, 1 + λ), 0)\n", + " H += diagm(vcat(-2λ, fill(-4λ, N - 3), -2λ), 1)\n", + " H += diagm(fill(λ, N-2), 2)\n", + " return H\n", + "end;" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "function hp_filter(y::AbstractVector{T}, λ::Real, mt::MatrixType, am::AssignMethod) where T <: Real\n", + " N = length(y)\n", + " H = create_matrix(T(λ), N, mt, am)\n", + " y_trend = H \\ y\n", + " y_cyclical = y - y_trend\n", + " return y_trend, y_cyclical\n", + "end;" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Apply HP filter" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*Note*: First run for each include compilation time, so don't take it seriously. Run a second time to get a better sense of the differences." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Looped dense\n", + " 5.267661 seconds (1.31 M allocations: 69.188 MiB, 6.57% gc time)\n", + " 0.048105 seconds (998 allocations: 954.137 KiB)\n", + "Looped sparse\n", + " 1.623298 seconds (250.93 k allocations: 13.624 MiB, 2.09% gc time)\n", + " 0.025615 seconds (1.21 k allocations: 255.492 KiB)\n", + "Vectorized dense\n", + " 1.128149 seconds (288.31 k allocations: 20.795 MiB, 3.79% gc time)\n", + " 0.024424 seconds (1.24 k allocations: 5.349 MiB)\n", + "Vectorized sparse\n", + " 3.065138 seconds (1.49 M allocations: 81.303 MiB, 3.46% gc time)\n", + " 0.014130 seconds (1.39 k allocations: 284.764 KiB)\n" + ] + } + ], + "source": [ + "data = readxlsheet(\"data/US_Data.xlsx\", \"Data\");\n", + "y = collect(data[4:end, 2]); # removing header\n", + "T = length(y);\n", + "\n", + "# Looped dense version\n", + "println(\"Looped dense\")\n", + "@time (ytr1600, yc1600) = hp_filter(Float64.(y), 1600, Dense(), Loop());\n", + "@time (ytr1e5, yc1e5) = hp_filter(Float64.(y), 1e5, Dense(), Loop());\n", + "\n", + "# Looped sparse version\n", + "println(\"Looped sparse\")\n", + "@time (ytr1600, yc1600) = hp_filter(Float64.(y), 1600, Sparse(), Loop());\n", + "@time (ytr1e5, yc1e5) = hp_filter(Float64.(y), 1e5, Sparse(), Loop());\n", + "\n", + "# Dense matrix version\n", + "println(\"Vectorized dense\")\n", + "@time (ytr1600, yc1600) = hp_filter(Float64.(y), 1600, Dense(), Vectorized());\n", + "@time (ytr1e5, yc1e5) = hp_filter(Float64.(y), 1e5, Dense(), Vectorized());\n", + "\n", + "# Sparse matrix version\n", + "println(\"Vectorized sparse\")\n", + "@time (ytr1600, yc1600) = hp_filter(Float64.(y), 1600, Sparse(), Vectorized());\n", + "@time (ytr1e5, yc1e5) = hp_filter(Float64.(y), 1e5, Sparse(), Vectorized());" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## With large data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Looped dense\n", + " 8.708596 seconds (27 allocations: 381.585 MiB, 5.26% gc time)\n", + "Looped sparse\n", + " " + ] + } + ], + "source": [ + "N_test = 5_000\n", + "y_test = randn(N_test)\n", + "\n", + "# Looped dense version\n", + "println(\"Looped dense\")\n", + "@time (ytrtest_del, yctest_del) = hp_filter(y_test, 1600, Dense(), Loop());\n", + "\n", + "# Looped sparse version\n", + "println(\"Looped sparse\")\n", + "@time (ytrtest_spl, yctest_spl) = hp_filter(y_test, 1600, Sparse(), Loop());\n", + "\n", + "# Dense matrix version\n", + "println(\"Vectorized dense\")\n", + "@time (ytrtest_dev, yctest_dev) = hp_filter(y_test, 1600, Dense(), Vectorized());\n", + "\n", + "# # Sparse matrix version\n", + "println(\"Vectorized sparse\")\n", + "@time (ytrtest_spv, yctest_spv) = hp_filter(y_test, 1600, Sparse(), Vectorized());" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot trend component\n", + "Now back to the actual assignment." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot([log.(ytr1600), log.(ytr1e5), log.(y)], lw = 2,\n", + " lab = [\"λ = 1600\" \"λ = 10⁵\" \"Raw Data\"],\n", + " title = \"log of real GDP\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot cyclical component" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "p = plot([yc1600, yc1e5], lw = 2,\n", + " lab = [\"λ = 1600\" \"λ = 10⁵\" \"Raw Data\"],\n", + " title = \"log of real GDP\")" + ] + } + ], + "metadata": { + "hide_input": false, + "kernelspec": { + "display_name": "Julia 0.6.0", + "language": "julia", + "name": "julia-0.6" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "0.6.1" + }, + "toc": { + "colors": { + "hover_highlight": "#DAA520", + "running_highlight": "#FF0000", + "selected_highlight": "#FFD700" + }, + "moveMenuLeft": true, + "nav_menu": { + "height": "161px", + "width": "253px" + }, + "navigate_menu": true, + "number_sections": true, + "sideBar": true, + "threshold": 4, + "toc_cell": false, + "toc_section_display": "block", + "toc_window_display": false + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}