diff --git a/questions_AD.ipynb b/questions_AD.ipynb new file mode 100644 index 0000000..6ca6b16 --- /dev/null +++ b/questions_AD.ipynb @@ -0,0 +1,270 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Integration exercises\n", + "\n", + "1. Fackler and Miranda exercise 5.1: Demand for a commodity is given by $q(p)=2p^{-0.5}$. The price of the good falls from 4 to 1. Compute the change in consumer surplus **analytically** using calculus and **numerically** using three different methods: a Gauss-Legendre rule, Monte-Carlo, and Quasi Monte Carlo.\n", + " 1. Start out by plotting the demand function. In that plot, label the axis properly ($p$ on the $x$ axis), and add 2 horizontal lines for the equilibrium quantities at both prices $p=1,p=4$. This should help you visualize the consumer surplus, and it should guide you in the analytical solution. The simplest way to solve this first part would be to edit this `IJulia` notebook, just adding a code cell below (for the plot), and then one for the analytic solution. You can get and see how to use `IJulia` [by clicking on this link](https://github.com/JuliaLang/IJulia.jl).\n", + "\n", + " To implement the integration rules, write 3 functions, one for each sub question, and detailed below. Let each function take as argument the number of integration points. Use $n=10,15,20$ points. Each of those functions should produce the result (of course), as well as a plot where we can see integration nodes vs function value. We want to get an idea where each method places the points, and how this might influence the different results. Ideally, the result would tell us how far the corresponding method is away from your analytic solution. *Ideally*, your plot would show us all the results at once, so we can easily compare across methods. Keyword *subplots*.\n", + " 1. Gauss-legendre rule. Note that you have to change the function domain to $[-1,1]$ first. This is achieved with the following transformation\n", + "$$ \\int_a^b f(x)\\,dx = \\frac{b-a}{2} \\int_{-1}^1 f\\left(\\frac{b-a}{2}x + \\frac{a+b}{2}\\right)dx $$\n", + " 1. monte-carlo, again taking $n$ as an argument.\n", + " 1. pseudo monte-carlo. Use a Sobol Sequence.\n", + "2. Fackler and Miranda exercise 5.5: A government stabilizes the supply of a commodity at $S=2$, but allows the price to be determined on the market. Domestic and export demand for the commodity are given by $D = \\tilde{\\theta_1} p^{-1}, X = \\tilde{\\theta_2} p^{-0.5}$, and where $\\log \\tilde{\\theta_1}$ and $\\log \\tilde{\\theta_1}$ are both normally distributed with means 0, variances 0.02 and 0.01 respectively, and covariance 0.01.\n", + " 1. compute the expected price $E[p]$ and expected variance $V[p]$ using a kronecker product rule with gauss-hermite grids (10) in each dimension of the shock.\n", + " 1. perform the same computation with a monte carlo integration scheme.\n", + " 1. Bonus question: same thing with quasi-monte carlo." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "#### Question 1. Part A." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "1\n", + "\n", + "\n", + "2\n", + "\n", + "\n", + "3\n", + "\n", + "\n", + "4\n", + "\n", + "\n", + "1.0\n", + "\n", + "\n", + "1.5\n", + "\n", + "\n", + "2.0\n", + "\n", + "\n", + "2.5\n", + "\n", + "\n", + "Demand curve\n", + "\n", + "\n", + "p\n", + "\n", + "\n", + "quantity\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "q(p)\n", + "\n", + "\n", + "\n", + "qty at p= 4\n", + "\n", + "\n" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Plot the demand function\n", + "using Plots\n", + "q(p) = 2*p^(-1/2)\n", + "\n", + "plot(q,0.5,5,label=\"q(p)\", title=\"Demand curve\", xaxis=\"p\", yaxis=\"quantity\")\n", + "hline!([1, 2], labels = [\"qty at p= 4\" \"qty at p'= 1\"], color = [\"red\", \"green\"])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Analytical solution" + ] + }, + { + "cell_type": "code", + "execution_count": 92, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Change in consumer surplus" + ] + }, + { + "data": { + "text/latex": [ + "$ \\Delta CS = \\int_1^{\\infty} q(p)\\, \\mathrm{d}p - \\int_4^{\\infty} q(p)\\, \\mathrm{d}p \n", + "= \\int_1^{4} q(p)\\,\\mathrm{d}p = \\int_1^{4} 2p^{-1/2}\\, \\mathrm{d}p\n", + "= \\int_1^{4} q(p)\\,\\mathrm{d}p = \\int_1^{4} 2p^{-1/2}\\, \\mathrm{d}p\n", + "= [4p^{1/2}]_1^4\n", + "= 4(4^{1/2}-1^{1/2})\n", + "= 4\n", + "$" + ], + "text/plain": [ + "L\"$ \\Delta CS = \\int_1^{\\infty} q(p)\\, \\mathrm{d}p - \\int_4^{\\infty} q(p)\\, \\mathrm{d}p \n", + "= \\int_1^{4} q(p)\\,\\mathrm{d}p = \\int_1^{4} 2p^{-1/2}\\, \\mathrm{d}p\n", + "= \\int_1^{4} q(p)\\,\\mathrm{d}p = \\int_1^{4} 2p^{-1/2}\\, \\mathrm{d}p\n", + "= [4p^{1/2}]_1^4\n", + "= 4(4^{1/2}-1^{1/2})\n", + "= 4\n", + "$\"" + ] + }, + "execution_count": 92, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using LaTeXStrings\n", + "\n", + "print(\"Change in consumer surplus\")\n", + "L\"\"\" \\Delta CS = \\int_1^{\\infty} q(p)\\, \\mathrm{d}p - \\int_4^{\\infty} q(p)\\, \\mathrm{d}p \n", + "= \\int_1^{4} q(p)\\,\\mathrm{d}p = \\int_1^{4} 2p^{-1/2}\\, \\mathrm{d}p\n", + "= \\int_1^{4} q(p)\\,\\mathrm{d}p = \\int_1^{4} 2p^{-1/2}\\, \\mathrm{d}p\n", + "= [4p^{1/2}]_1^4\n", + "= 4(4^{1/2}-1^{1/2})\n", + "= 4\n", + "\"\"\"\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 0.6.2", + "language": "julia", + "name": "julia-0.6" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "0.6.2" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/src/HWintegration.jl b/src/HWintegration.jl index 68ef546..e163b4a 100644 --- a/src/HWintegration.jl +++ b/src/HWintegration.jl @@ -1,7 +1,7 @@ module HWintegration - const A_SOL = -18 # enter your analytic solution. -18 is wrong. + const A_SOL = 4 # imports using FastGaussQuadrature @@ -17,51 +17,83 @@ module HWintegration srand(12345) # demand function - - # gauss-legendre adjustment factors for map change - - # eqm condition for question 2 + q(p) = 2p^(-1/2) # makes a plot for questions 1 function plot_q1() - + plot(q,0.5,5,label="q(p)", title="Demand curve", xaxis="p", yaxis="quantity") + hline!([1, 2], label=["qty at p= 4", "qty at p'= 1"], color = ["red", "green"]) end - function question_1b(n) - + #compute node points and weights + x = gausslegendre(n) + #Change function interval + nodes_GL = 3/2 * x[1] + 5/2 + values_GL = [q(i) for i in nodes_GL] + #Compute approximate integral + GL_int = 3/2 * sum([x[2][i] * values_GL[i] for i in 1:n]) + deviation = A_SOL - GL_int + # Plot graph + global GL_graph = scatter(nodes_GL, values_GL, title="Gauss Legendre, n=$n") + println("Gauss-Legendre method with $n nodes yields a change of $GL_int. + Deviation from analytical solution is $deviation") end - function question_1c(n) + nodes_MC = rand(Uniform(1,4), n) + values_MC = [q(x) for x in nodes_MC] + MC_int = 3/n * sum(values_MC) + deviation = A_SOL - MC_int + + #graph + global MC_graph = scatter(nodes_MC, values_MC, title="Monte Carlo, n=$n") - + println("Monte Carlo method with $n nodes yields a change of $MC_int. + Deviation from analytical solution is $deviation") end function question_1d(n) - - - + # Generate n points in one dimension (rescaled to be evenly spread on [1,4] interval) + s = SobolSeq(1) + nodes_qMC = 3*[[next(s) for i in 1:n][x][1] for x in 1:n] + 1 + values_qMC = [q(i) for i in nodes_qMC] + qMC_int = 3/n * sum(values_qMC) + deviation = A_SOL - qMC_int + + #graph + global qMC_graph = scatter(nodes_qMC, values_qMC, title="Quasi Monte Carlo, n=$n") + + println("Quasi Monte Carlo method with $n nodes yields a change of $qMC_int. + Deviation from analytical solution is $deviation") end # question 2 + # eqm condition for question 2 + # theta_1*p^(-1) + theta_2*p^(-1/2) = 2 + eqm(theta, supply, p) = theta[1]*p + theta[2]*p - supply + function question_2a(n) + mu = [0, 0] + sigma = [0.02 0.01; 0.01 0.01] - + # Gauss-Hermite nodes + #nodes = Any[] + #push!(nodes,repeat(gausshermite(n),inner=[1],outer=[n])) # dim1 + #push!(nodes,repeat(rules["hermite"][1],inner=[3],outer=[3])) # dim2 + #weights = kron(rules["hermite"][2],kron(rules["hermite"][2],rules["hermite"][2])) + #df = hcat(DataFrame(weights=weights),DataFrame(nodes,[:dim1,:dim2])) + println("I could not solve this question") end function question_2b(n) - - + println("I could not solve this question") end function question_2bonus(n) - - - - + println("I could not solve this question") end # function to run all questions @@ -72,13 +104,13 @@ module HWintegration for n in (10,15,20) info("============================") info("now showing results for n=$n") - # info("question 1b:") - # question_1b(n) # make sure your function prints some kind of result! - # info("question 1c:") - # question_1c(n) - # info("question 1d:") - # question_1d(n) - # println("") + info("question 1b:") + question_1b(n) # make sure your function prints some kind of result! + info("question 1c:") + question_1c(n) + info("question 1d:") + question_1d(n) + println("") info("question 2a:") q2 = question_2a(n) println(q2)