diff --git a/.ipynb_checkpoints/questions-checkpoint.ipynb b/.ipynb_checkpoints/questions-checkpoint.ipynb
new file mode 100644
index 0000000..4fe87c5
--- /dev/null
+++ b/.ipynb_checkpoints/questions-checkpoint.ipynb
@@ -0,0 +1,233 @@
+{
+ "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": {},
+ "source": [
+ "## (A) Plot"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n"
+ ]
+ },
+ "execution_count": 1,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "using Plots\n",
+ "gr()\n",
+ "x=0.5:10\n",
+ "plot(x, 2 ./ sqrt.(x))\n",
+ "vline!([4], color=:green,lw=1,label=\"p = 4\")\n",
+ "vline!([1], color=:orange,lw=1,label=\"p = 1\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## (A) Analytical solution\n",
+ "\n",
+ "$$q(p) = 2p^{0.5}$$\n",
+ "\n",
+ "The change in consumer surplus is the difference between the consumer surplus after and before the change in prices. Let denote $\\Delta$ that number. \n",
+ "\n",
+ "$\\begin{array}{lcl}\n",
+ "\\Delta & = & \\displaystyle\\int_1^\\infty q(p)dp - \\displaystyle\\int_4^\\infty q(p)dp \\\\\n",
+ "\\Delta & = & \\displaystyle\\int_1^4 q(p)dp \\\\\n",
+ "\\Delta & = & \\left[4\\sqrt{p}\\right]_1^4 \\\\\n",
+ "\\Delta & = & 4 \\\\\n",
+ "\\end{array}$\n",
+ "\n",
+ "The change in consumer surplus is $4$. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "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/.travis.yml b/.travis.yml
index 7a3fe23..d404cca 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,10 +1,13 @@
language: julia
-os:
- - linux
julia:
- 0.6
notifications:
email: false
+testscript:
+ julia -e 'Pkg.build("HWintegration")'
+if [ -f test/runtests.jl ]; then
+ julia --check-bounds=yes -e 'Pkg.test("HWintegration", coverage=true)'
+fi
#script: # the default script is equivalent to the following
# - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
# - julia -e 'Pkg.clone(pwd()); Pkg.build("Example"); Pkg.test("Example"; coverage=true)';
diff --git a/questions.ipynb b/questions.ipynb
index bead3f5..4fe87c5 100644
--- a/questions.ipynb
+++ b/questions.ipynb
@@ -20,19 +20,204 @@
" 1. Bonus question: same thing with quasi-monte carlo."
]
},
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## (A) Plot"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n"
+ ]
+ },
+ "execution_count": 1,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "using Plots\n",
+ "gr()\n",
+ "x=0.5:10\n",
+ "plot(x, 2 ./ sqrt.(x))\n",
+ "vline!([4], color=:green,lw=1,label=\"p = 4\")\n",
+ "vline!([1], color=:orange,lw=1,label=\"p = 1\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## (A) Analytical solution\n",
+ "\n",
+ "$$q(p) = 2p^{0.5}$$\n",
+ "\n",
+ "The change in consumer surplus is the difference between the consumer surplus after and before the change in prices. Let denote $\\Delta$ that number. \n",
+ "\n",
+ "$\\begin{array}{lcl}\n",
+ "\\Delta & = & \\displaystyle\\int_1^\\infty q(p)dp - \\displaystyle\\int_4^\\infty q(p)dp \\\\\n",
+ "\\Delta & = & \\displaystyle\\int_1^4 q(p)dp \\\\\n",
+ "\\Delta & = & \\left[4\\sqrt{p}\\right]_1^4 \\\\\n",
+ "\\Delta & = & 4 \\\\\n",
+ "\\end{array}$\n",
+ "\n",
+ "The change in consumer surplus is $4$. "
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "collapsed": true
- },
+ "metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
- "display_name": "Julia 0.6.0",
+ "display_name": "Julia 0.6.2",
"language": "julia",
"name": "julia-0.6"
},
@@ -40,7 +225,7 @@
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
- "version": "0.6.0"
+ "version": "0.6.2"
}
},
"nbformat": 4,
diff --git a/src/HWintegration.jl b/src/HWintegration.jl
index 68ef546..fadab83 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
@@ -10,6 +10,7 @@ module HWintegration
using Plots
using Distributions
using NullableArrays
+ using DataFrames
using DataStructures # OrderedDict
@@ -18,50 +19,202 @@ module HWintegration
# demand function
+ q(x)=2/sqrt(x)
+
# gauss-legendre adjustment factors for map change
+ a = 4
+ b = 1
+ coef1 = (a-b)/2
+ coef2 = (a+b)/2
+
# eqm condition for question 2
+ # S = X + D
+ # theta_1/x + theta_2/sqrt(x) - 2 = 0
+
# makes a plot for questions 1
function plot_q1()
-
+ gr()
+ x=0.5:10
+ plot(x, 2 ./ sqrt.(x))
+ vline!([4], color=:green,lw=1,label="p = 4")
+ vline!([1], color=:orange,lw=1,label="p = 1")
end
- function question_1b(n)
+ function question_1b(n)
+ q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2)
+ rule = gausslegendre(n)
+ nodes = values(rule)[1,]
+ weights = values(rule)[2,]
+ weighted_q = weights .* map(q_trans, nodes)
+ approx = sum(weighted_q)
+ println("When n = $n :")
+ println("The approximated change in consumer surplus is $approx.")
+ println("The distance between the true result and the approximation is ", A_SOL - approx)
+ scatter(nodes, map(q_trans, nodes), label = "Gauss-Legendre", xlab ="Integration nodes", ylab="Function Value" )
+ plot!(q_trans, label = "Scaled Demand Function")
+ end
- end
+ function question_1c(n)
- function question_1c(n)
+ q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2)
+ x = sort(rand(n))
+ y = map(q_trans, x)
+ approx = sum(y)*(1+1)/n
+ println("When n = $n :")
+ println("The approximated change in consumer surplus is ", approx)
+ println("The distance between the true result and the approximation is ", A_SOL - approx)
+ scatter!(x, y, label = "Monte Carlo", xlab ="Integration nodes", ylab="Function Value" )
-
- end
+ end
- function question_1d(n)
-
+ function question_1d(n)
- end
+ #= We define a function Sobol that generates the sobol
+ sequence in a usable way =#
- # question 2
+ function sobol(j)
+ s = ScaledSobolSeq(1, [-1.0], [1.0])
+ seq = hcat([next(s) for i in 1:j]...)
+ a = []
+ for i in 1:j
+ push!(a, seq[1, i])
+ end
+ return a
+ end
- function question_2a(n)
+ q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2)
+
+ rand = sobol(n)
+ x = sort(rand)
+ y = map(q_trans, x)
+
+ approx = sum(y)*(1+1)/n
+ println("When n = $n :")
+ println("The approximated change in consumer surplus is $approx.")
+ println("The distance between the true result and the approximation is ", A_SOL - approx)
+ scatter!(x, y, label = "Q. Monte Carlo", xlab ="Integration nodes", ylab="Function Value" )
+
+ end
-
+ # question 2
+ function question_2a(n)
+ rule = gausshermite(n)
+ nodes = Any[]
+ push!(nodes,repeat(rule[1],inner=[1],outer=[n]))
+ push!(nodes,repeat(rule[1],inner=[n],outer=[1]))
+ #Making the nodes being of variance 1
+ nodes[1] = nodes[1]/sqrt(var(nodes[1]))
+ nodes[2] = nodes[2]/sqrt(var(nodes[2]))
+ weights = kron(rule[2],rule[2])
+ points = DataFrame(weights = weights)
+ points[:dim1] = nodes[1]
+ points[:dim2] = nodes[2]
+ theta_unnorm = hcat(points[2], points[3])
+
+ #= Making log(theta1) and log(theta2) such that
+ var(log(theta1)) = 0.01
+ var(log(theta2)) = 0.02
+ cov(log(theta1), log(theta2)) = 0.01 =#
+ dim = n^2
+ SIG = [sqrt(0.01) 0 ; sqrt(0.01) sqrt(0.01)]
+ theta = (SIG*theta_unnorm')'
+ points[:theta_1] = exp.(theta[1:dim,])
+ points[:theta_2] = exp.(theta[dim+1:dim*2,])
+
+ # Generating price vector
+
+ price = []
+ for i in 1:dim
+ a = points[:theta_1][i]
+ b = points[:theta_2][i]
+ push!(price, fzero(x -> a/x + b/sqrt(x) - 2, 0, 10))
+ end
+
+ points[:price] = price
+
+ approx_E = sum(points[:price] .* points[:weights])
+ approx_V = sum((points[:price] .^ 2) .* points[:weights])
+ println("When n = $n :")
+ println("The approximated expected price is $approx_E.")
+ println("The approximated expected variance is $approx_V")
end
function question_2b(n)
+ L = 0.1 * [1/sqrt(2) 0; 1/sqrt(2) sqrt(2)]
+ theta_1 = randn(n)
+ theta_2 = randn(n)
+ theta_unnorm = hcat(theta_1, theta_2)
+ theta_norm = (L*theta_unnorm')'
+ theta_1_norm = [exp.(theta_norm[1:n,])]
+ theta_2_norm = [exp.(theta_norm[n+1:n*2,])]
+ points = DataFrame(theta_1_norm = theta_1_norm[1], theta_2_norm = theta_2_norm[1])
+
+ price = []
+ for i in 1:n
+ a = points[:theta_1_norm][i]
+ b = points[:theta_2_norm][i]
+ if (a<0) && (b<0)
+ push!(price,0)
+ else
+ push!(price, fzero(x -> a/x + b/sqrt(x) - 2, 0, 10))
+ end
+ end
+ points[:price] = price
+ approx_E = sum(points[:price])/n
+ approx_V = sum( (points[:price]- approx_E) .^ 2 ) / n
+ println("When n = $n :")
+ println("The approximated expected price is $approx_E.")
+ println("The approximated expected variance is $approx_V")
-
end
function question_2bonus(n)
-
-
-
+ function sobol(j)
+ s = ScaledSobolSeq(1, [-1.0], [1.0])
+ seq = hcat([next(s) for i in 1:j]...)
+ a = []
+ for i in 1:j
+ push!(a, seq[1, i])
+ end
+ return a
+ end
+
+ nodes = sobol(n)
+ theta_unnorm = []
+ push!(theta_unnorm,repeat(nodes,inner=[1],outer=[n]))
+ push!(theta_unnorm,repeat(nodes,inner=[n],outer=[1]))
+ L = 0.1 * [1/sqrt(2) 0; 1/sqrt(2) sqrt(2)]
+ points = DataFrame(theta_1 = theta_unnorm[1], theta_2 = theta_unnorm[2])
+ theta = hcat(points[1], points[2])
+ theta_norm = (L*theta')'
+ dim = n^2
+ theta_1_norm = [theta_norm[1:dim,]]
+ theta_2_norm = [theta_norm[dim+1:dim*2,]]
+ points[:theta_1_norm] = exp.(theta_1_norm[1])
+ points[:theta_2_norm] = exp.(theta_2_norm[1])
+ price = []
+ for i in 1:dim
+ a = points[:theta_1_norm][i]
+ b = points[:theta_2_norm][i]
+ if (a<0) && (b<0)
+ push!(price,0)
+ else
+ push!(price, fzero(x -> a/x + b/sqrt(x) - 2, 0, 10))
+ end
+ end
+ points[:price] = price
+ approx_E = sum(points[:price])/dim
+ approx_V = sum( (points[:price]- approx_E) .^ 2 ) / dim
+ println("When n = $n :")
+ println("The approximated expected price is $approx_E.")
+ println("The approximated expected variance is $approx_V")
end
# function to run all questions
@@ -69,26 +222,31 @@ module HWintegration
info("Running all of HWintegration")
info("question 1:")
plot_q1()
- 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 2a:")
+ info("question 2:")
+ plotlyjs()
+ plots = []
+ for n in (10,15,20)
+ info("============================")
+ info("now showing results for n=$n")
+ info("question 1b (Gauss Legendre):")
+ question_1b(n)
+ info("question 1c (Monte Carlo):")
+ question_1c(n)
+ info("question 1d (Quasi Monte Carlo):")
+ push!(plots, question_1d(n))
+
+ info("question 2:")
+ info("question 2a (Gauss Hermite)")
q2 = question_2a(n)
println(q2)
- info("question 2b:")
+ info("question 2b (Monte Carlo):")
q2b = question_2b(n)
println(q2b)
- info("bonus question: Quasi monte carlo:")
+ info("bonus question (Quasi monte carlo):")
q2bo = question_2bonus(n)
println(q2bo)
end
+ plot(plots[1], plots[2], plots[3])
end
info("end of HWintegration")
diff --git a/test/runtests.jl b/test/runtests.jl
index 1bb9cb0..362e38b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -10,14 +10,17 @@ using Base.Test
@testset "check gauss adjustment" begin
+ @test HWintegration.coef1 == 1.5
+ @test HWintegration.coef2 == 2.5
end
@testset "eqm condition for Q2" begin
+ @test 1==1
end
end
@testset "Testing Results" begin
-
+ @test 1 ==1
end
-
+
end