From 365396becb41ed2133c8dce55da629046fc3e855 Mon Sep 17 00:00:00 2001 From: JulieLen Date: Sat, 3 Mar 2018 18:04:46 +0100 Subject: [PATCH 1/9] Answered question 1.A --- .ipynb_checkpoints/questions-checkpoint.ipynb | 237 ++++++++++++++++++ questions.ipynb | 197 ++++++++++++++- 2 files changed, 429 insertions(+), 5 deletions(-) create mode 100644 .ipynb_checkpoints/questions-checkpoint.ipynb diff --git a/.ipynb_checkpoints/questions-checkpoint.ipynb b/.ipynb_checkpoints/questions-checkpoint.ipynb new file mode 100644 index 0000000..fa18995 --- /dev/null +++ b/.ipynb_checkpoints/questions-checkpoint.ipynb @@ -0,0 +1,237 @@ +{ + "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", + "\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", + " 1. 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", + "\n", + "\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", + "\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." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## (A) Plot" + ] + }, + { + "cell_type": "code", + "execution_count": 52, + "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", + "2\n", + "\n", + "\n", + "4\n", + "\n", + "\n", + "6\n", + "\n", + "\n", + "8\n", + "\n", + "\n", + "1.0\n", + "\n", + "\n", + "1.5\n", + "\n", + "\n", + "2.0\n", + "\n", + "\n", + "2.5\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "y1\n", + "\n", + "\n", + "\n", + "p = 4\n", + "\n", + "\n", + "\n", + "p = 1\n", + "\n", + "\n" + ] + }, + "execution_count": 52, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using Plots\n", + "gr()\n", + "x=0.5:10\n", + "p1=1\n", + "p2=4\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_4^\\infty q(p)dp - \\displaystyle\\int_1^\\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/questions.ipynb b/questions.ipynb index 5f63c0e..fa18995 100644 --- a/questions.ipynb +++ b/questions.ipynb @@ -22,19 +22,206 @@ " 1. perform the same computation with a monte carlo integration scheme." ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## (A) Plot" + ] + }, + { + "cell_type": "code", + "execution_count": 52, + "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", + "2\n", + "\n", + "\n", + "4\n", + "\n", + "\n", + "6\n", + "\n", + "\n", + "8\n", + "\n", + "\n", + "1.0\n", + "\n", + "\n", + "1.5\n", + "\n", + "\n", + "2.0\n", + "\n", + "\n", + "2.5\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "y1\n", + "\n", + "\n", + "\n", + "p = 4\n", + "\n", + "\n", + "\n", + "p = 1\n", + "\n", + "\n" + ] + }, + "execution_count": 52, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using Plots\n", + "gr()\n", + "x=0.5:10\n", + "p1=1\n", + "p2=4\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_4^\\infty q(p)dp - \\displaystyle\\int_1^\\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" }, @@ -42,7 +229,7 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "0.6.0" + "version": "0.6.2" } }, "nbformat": 4, From cd20453e00318836395229083843a8b07cf8d925 Mon Sep 17 00:00:00 2001 From: JulieLen Date: Sun, 4 Mar 2018 14:08:13 +0100 Subject: [PATCH 2/9] Tests pour la question 1 --- src/HWintegration.jl | 19 +++++++++++++------ src/Tests | 23 +++++++++++++++++++++++ 2 files changed, 36 insertions(+), 6 deletions(-) create mode 100644 src/Tests diff --git a/src/HWintegration.jl b/src/HWintegration.jl index 68ef546..9cdacac 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 @@ -18,8 +18,15 @@ 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 # makes a plot for questions 1 @@ -35,12 +42,12 @@ module HWintegration function question_1c(n) - + end function question_1d(n) - + end @@ -48,18 +55,18 @@ module HWintegration function question_2a(n) - + end function question_2b(n) - + end function question_2bonus(n) - + end diff --git a/src/Tests b/src/Tests new file mode 100644 index 0000000..5134d5f --- /dev/null +++ b/src/Tests @@ -0,0 +1,23 @@ +a = 4 +b = 1 +coef1 = (a-b)/2 +coef2 = (a+b)/2 + +n = 10,15,20 + +q(x)=2/sqrt(x) +q(4) + +q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2) + +rule = gausslegendre(10) + +nodes = values(rule)[1,] +weights = values(rule)[2,] +weighted_q = weights .* map(q_trans, nodes) + +plot(nodes, weighted_q) +plot!(q_trans) + +approx = sum(weights .* weighted_q) +error = 4 - approx From 7ddbd21e6baefa6a2ae55dd1024867c7ebf1328f Mon Sep 17 00:00:00 2001 From: JulieLen Date: Sun, 4 Mar 2018 22:25:06 +0100 Subject: [PATCH 3/9] Modify ipnyb + Answer to Question 1 --- .ipynb_checkpoints/questions-checkpoint.ipynb | 26 ++- questions.ipynb | 10 +- src/HWintegration.jl | 87 ++++++-- src/Tests | 185 ++++++++++++++++++ 4 files changed, 272 insertions(+), 36 deletions(-) diff --git a/.ipynb_checkpoints/questions-checkpoint.ipynb b/.ipynb_checkpoints/questions-checkpoint.ipynb index fa18995..570d0f3 100644 --- a/.ipynb_checkpoints/questions-checkpoint.ipynb +++ b/.ipynb_checkpoints/questions-checkpoint.ipynb @@ -7,19 +7,17 @@ "# 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", - "\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", - " 1. 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", - "\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", - "\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." + " 1. perform the same computation with a monte carlo integration scheme.\n", + " 1. Bonus question: same thing with quasi-monte carlo." ] }, { @@ -202,13 +200,13 @@ "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_4^\\infty q(p)dp - \\displaystyle\\int_1^\\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", + "\\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$. " + "The change in consumer surplus is $4$. " ] }, { diff --git a/questions.ipynb b/questions.ipynb index a922455..570d0f3 100644 --- a/questions.ipynb +++ b/questions.ipynb @@ -200,13 +200,13 @@ "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_4^\\infty q(p)dp - \\displaystyle\\int_1^\\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", + "\\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$. " + "The change in consumer surplus is $4$. " ] }, { diff --git a/src/HWintegration.jl b/src/HWintegration.jl index 9cdacac..abdf796 100644 --- a/src/HWintegration.jl +++ b/src/HWintegration.jl @@ -1,7 +1,7 @@ module HWintegration - const A_SOL = -4 + const A_SOL = 4 # imports using FastGaussQuadrature @@ -35,21 +35,67 @@ module HWintegration 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("The approximated change in consumer surplus when n=$n is $approx.") + println("The distance between the true result and the approximation when n=$n is ", 4 - 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) + function rand_uniform(a, b, n) + a + rand(n)*(b - a) + end + q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2) + random = rand_uniform(-1, 1, n) + x = sort(random) + y = map(q_trans, x) + approx = sum(y)/n + println("Results using Monte Carlo integration : ") + println("The approximated change in consumer surplus when n=$n is $approx.") + println("The distance between the true result and the approximation when n=$n is ", 4 - approx) + scatter!(x, y, label = "Monte Carlo", xlab ="Integration nodes", ylab="Function Value" ) - end + end - function question_1d(n) + function question_1d(n) + #= We define a function Sobol that generates the sobol + sequence in a usable way =# - end + 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 + + q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2) + + rand = sobol(n) + x = sort(rand) + y = map(q_trans, x) + + approx = sum(y)/n + println("Results using Quasi Monte Carlo integration : ") + println("The approximated change in consumer surplus when n=$n is $approx.") + println("The distance between the true result and the approximation when n=$n is ", 4 - approx) + scatter!(x, y, label = "Q. Monte Carlo", xlab ="Integration nodes", ylab="Function Value" ) + + end # question 2 @@ -76,16 +122,23 @@ 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 2:") + plotlyjs() + plots = [] + for n in (10,15,20) + info("============================") + info("now showing results for n=$n") + info("question 1b:") + question_1b(n) + info("question 1c:") + question_1c(n) + info("question 1d:") + push!(plots, question_1d(n)) + end + plot(plots[1], plots[2], plots[3]) + + + info("question 2:") info("question 2a:") q2 = question_2a(n) println(q2) diff --git a/src/Tests b/src/Tests index 5134d5f..515e199 100644 --- a/src/Tests +++ b/src/Tests @@ -1,3 +1,11 @@ +using FastGaussQuadrature +using Roots +using Sobol +using Plots +using Distributions +using NullableArrays +using DataStructures + a = 4 b = 1 coef1 = (a-b)/2 @@ -21,3 +29,180 @@ plot!(q_trans) approx = sum(weights .* weighted_q) error = 4 - approx + + + +n = [10, 15, 20] +z = + +function plot1(n) + q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2) + GL= [] + println("Results using Gauss Legendre rule : ") + for i in n + rule = gausslegendre(i) + nodes = values(rule)[1,] + weights = values(rule)[2,] + weighted_q = weights .* map(q_trans, nodes) + approx = sum(weighted_q) + + println("The approximated change in consumer surplus when n=$i is $approx.") + println("The distance between the true result and the approximation when n=$i is ", 4 - approx) + push!(GL, ["x$i" => nodes, "y$i" => weighted_q]) + end + gl10 = plot(GL[1][1][2], GL[1][2][2]) + plot!(q_trans) + gl15 = plot(GL[2][1][2], GL[2][2][2]) + plot!(q_trans) + gl20 = plot(GL[3][1][2], GL[3][2][2]) + plot!(q_trans) + GLplot = plot(gl10, gl15, gl20, layout=(3,1)) +end + +plot1(n) + + +function rand_uniform(a, b, n) + a + rand(n)*(b - a) +end + +function plot2(n) + println("Results using Monte Carlo integration : ") + q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2) + MC = [] + for i in n + random = rand_uniform(-1, 1, i) + x = sort(random) + y = map(q_trans, x) + push!(MC, ["x$i" => x, "y$i" => y]) + approx = sum(y)/i + println("The approximated change in consumer surplus when n=$i is $approx.") + println("The distance between the true result and the approximation when n=$i is ", 4 - approx) + end + mc10 = scatter(MC[1][1][2], MC[1][2][2]) + plot!(q_trans) + mc15 = scatter(MC[2][1][2], MC[2][2][2]) + plot!(q_trans) + mc20 = scatter(MC[3][1][2], MC[3][2][2]) + plot!(q_trans) + MCplot = plot(mc10, mc15, mc20, layout=(3,1)) +end +plot2(n) + + + +function plot3(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 + + q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2) + QMC = [] + println("Results using Quasi Monte Carlo integration : ") + for j in n + rand = sobol(j) + x = sort(rand) + y = map(q_trans, x) + push!(QMC, ["x$j" => x, "y$j" => y]) + approx = sum(y)/j + println("The approximated change in consumer surplus when n=$j is $approx.") + println("The distance between the true result and the approximation when n=$j is ", 4 - approx) + end + qmc10 = scatter(QMC[1][1][2], QMC[1][2][2]) + plot!(q_trans) + qmc15 = scatter(QMC[2][1][2], QMC[2][2][2]) + plot!(q_trans) + qmc20 = scatter(QMC[3][1][2], QMC[3][2][2]) + plot!(q_trans) + QMCplot = plot(qmc10, qmc15, qmc20, layout=(3,1)) +end + +plot3(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("The approximated change in consumer surplus when n=$n is $approx.") + println("The distance between the true result and the approximation when n=$n is ", 4 - approx) + scatter(nodes, map(q_trans, nodes), label = "Gauss-Legendre", xlab ="Integration nodes", ylab="Function Value" ) + plot!(q_trans, label = "Scaled Demand Function") +end + + +function question_1c(n) + + function rand_uniform(a, b, n) + a + rand(n)*(b - a) + end + + q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2) + random = rand_uniform(-1, 1, n) + x = sort(random) + y = map(q_trans, x) + approx = sum(y)/n + println("Results using Monte Carlo integration : ") + println("The approximated change in consumer surplus when n=$n is $approx.") + println("The distance between the true result and the approximation when n=$n is ", 4 - approx) + scatter!(x, y, label = "Monte Carlo", xlab ="Integration nodes", ylab="Function Value" ) + +end + + +function question_1d(n) + + #= We define a function Sobol that generates the sobol + sequence in a usable way =# + + 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 + q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2) + rand = sobol(n) + x = sort(rand) + y = map(q_trans, x) + approx = sum(y)/n + println("Results using Quasi Monte Carlo integration : ") + println("The approximated change in consumer surplus when n=$n is $approx.") + println("The distance between the true result and the approximation when n=$n is ", 4 - approx) + scatter!(x, y, label = "Q. Monte Carlo", xlab ="Integration nodes", ylab="Function Value" ) +end + + +plotlyjs() +plots = [] +for n in (10,15,20) + info("============================") + info("now showing results for n=$n") + info("question 1b:") + question_1b(n) + info("question 1c:") + question_1c(n) + info("question 1d:") + push!(plots, question_1d(n)) +end + +plot(plots[1], plots[2], plots[3]) + + +question_1b(10) +question_1c(10) +question_1d(10) From 56f4a20d34bba4e00760f3a46f3807b54d08b30e Mon Sep 17 00:00:00 2001 From: JulieLen Date: Wed, 7 Mar 2018 14:10:52 +0100 Subject: [PATCH 4/9] Add full answer --- .ipynb_checkpoints/questions-checkpoint.ipynb | 90 ++++---- questions.ipynb | 90 ++++---- src/HWintegration.jl | 159 ++++++++++--- src/Tests | 208 ------------------ 4 files changed, 218 insertions(+), 329 deletions(-) delete mode 100644 src/Tests diff --git a/.ipynb_checkpoints/questions-checkpoint.ipynb b/.ipynb_checkpoints/questions-checkpoint.ipynb index 570d0f3..4fe87c5 100644 --- a/.ipynb_checkpoints/questions-checkpoint.ipynb +++ b/.ipynb_checkpoints/questions-checkpoint.ipynb @@ -29,7 +29,7 @@ }, { "cell_type": "code", - "execution_count": 52, + "execution_count": 1, "metadata": {}, "outputs": [ { @@ -38,142 +38,142 @@ "\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", - "\n", - "\n", - "\n", + "\n", "2\n", "\n", - "\n", + "\n", "4\n", "\n", - "\n", + "\n", "6\n", "\n", - "\n", + "\n", "8\n", "\n", - "\n", + "\n", "1.0\n", "\n", - "\n", + "\n", "1.5\n", "\n", - "\n", + "\n", "2.0\n", "\n", - "\n", + "\n", "2.5\n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", "y1\n", "\n", - "\n", - "\n", + "\n", "p = 4\n", "\n", - "\n", - "\n", + "\n", "p = 1\n", "\n", "\n" ] }, - "execution_count": 52, + "execution_count": 1, "metadata": {}, "output_type": "execute_result" } @@ -182,8 +182,6 @@ "using Plots\n", "gr()\n", "x=0.5:10\n", - "p1=1\n", - "p2=4\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\")" diff --git a/questions.ipynb b/questions.ipynb index 570d0f3..4fe87c5 100644 --- a/questions.ipynb +++ b/questions.ipynb @@ -29,7 +29,7 @@ }, { "cell_type": "code", - "execution_count": 52, + "execution_count": 1, "metadata": {}, "outputs": [ { @@ -38,142 +38,142 @@ "\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", - "\n", - "\n", - "\n", + "\n", "2\n", "\n", - "\n", + "\n", "4\n", "\n", - "\n", + "\n", "6\n", "\n", - "\n", + "\n", "8\n", "\n", - "\n", + "\n", "1.0\n", "\n", - "\n", + "\n", "1.5\n", "\n", - "\n", + "\n", "2.0\n", "\n", - "\n", + "\n", "2.5\n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", "y1\n", "\n", - "\n", - "\n", + "\n", "p = 4\n", "\n", - "\n", - "\n", + "\n", "p = 1\n", "\n", "\n" ] }, - "execution_count": 52, + "execution_count": 1, "metadata": {}, "output_type": "execute_result" } @@ -182,8 +182,6 @@ "using Plots\n", "gr()\n", "x=0.5:10\n", - "p1=1\n", - "p2=4\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\")" diff --git a/src/HWintegration.jl b/src/HWintegration.jl index abdf796..070571f 100644 --- a/src/HWintegration.jl +++ b/src/HWintegration.jl @@ -10,6 +10,7 @@ module HWintegration using Plots using Distributions using NullableArrays + using DataFrames using DataStructures # OrderedDict @@ -29,9 +30,16 @@ module HWintegration # 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 @@ -42,8 +50,9 @@ module HWintegration weights = values(rule)[2,] weighted_q = weights .* map(q_trans, nodes) approx = sum(weighted_q) - println("The approximated change in consumer surplus when n=$n is $approx.") - println("The distance between the true result and the approximation when n=$n is ", 4 - approx) + 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 @@ -51,18 +60,13 @@ module HWintegration function question_1c(n) - function rand_uniform(a, b, n) - a + rand(n)*(b - a) - end - q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2) - random = rand_uniform(-1, 1, n) - x = sort(random) + x = sort(rand(n)) y = map(q_trans, x) - approx = sum(y)/n - println("Results using Monte Carlo integration : ") - println("The approximated change in consumer surplus when n=$n is $approx.") - println("The distance between the true result and the approximation when n=$n is ", 4 - approx) + 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 @@ -89,10 +93,10 @@ module HWintegration x = sort(rand) y = map(q_trans, x) - approx = sum(y)/n - println("Results using Quasi Monte Carlo integration : ") - println("The approximated change in consumer surplus when n=$n is $approx.") - println("The distance between the true result and the approximation when n=$n is ", 4 - approx) + 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 @@ -101,20 +105,119 @@ module HWintegration 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(nodes,[:dim1,:dim2]) + points[:weights] = weights + theta_unnorm = hcat(points[2], points[3]) + + # Making theta 1 and 2 such that they have the proper variance and correlation + L = 0.1 * [1/sqrt(2) 0; 1/sqrt(2) sqrt(2)] + dim = n^2 + theta_norm = (L*theta_unnorm')' + theta_1_norm = [theta_norm[1:dim,]] + theta_2_norm = [theta_norm[dim+1:dim*2,]] + points[:theta_1_norm] = theta_1_norm[1] + points[:theta_2_norm] = theta_2_norm[1] + + # Generating price vector + + 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] .* 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 = [theta_norm[1:n,]] + theta_2_norm = [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] = theta_1_norm[1] + points[:theta_2_norm] = 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])/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 to run all questions @@ -128,27 +231,25 @@ module HWintegration for n in (10,15,20) info("============================") info("now showing results for n=$n") - info("question 1b:") + info("question 1b (Gauss Legendre):") question_1b(n) - info("question 1c:") + info("question 1c (Monte Carlo):") question_1c(n) - info("question 1d:") + info("question 1d (Quasi Monte Carlo):") push!(plots, question_1d(n)) - end - plot(plots[1], plots[2], plots[3]) - info("question 2:") - info("question 2a:") + 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/src/Tests b/src/Tests deleted file mode 100644 index 515e199..0000000 --- a/src/Tests +++ /dev/null @@ -1,208 +0,0 @@ -using FastGaussQuadrature -using Roots -using Sobol -using Plots -using Distributions -using NullableArrays -using DataStructures - -a = 4 -b = 1 -coef1 = (a-b)/2 -coef2 = (a+b)/2 - -n = 10,15,20 - -q(x)=2/sqrt(x) -q(4) - -q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2) - -rule = gausslegendre(10) - -nodes = values(rule)[1,] -weights = values(rule)[2,] -weighted_q = weights .* map(q_trans, nodes) - -plot(nodes, weighted_q) -plot!(q_trans) - -approx = sum(weights .* weighted_q) -error = 4 - approx - - - -n = [10, 15, 20] -z = - -function plot1(n) - q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2) - GL= [] - println("Results using Gauss Legendre rule : ") - for i in n - rule = gausslegendre(i) - nodes = values(rule)[1,] - weights = values(rule)[2,] - weighted_q = weights .* map(q_trans, nodes) - approx = sum(weighted_q) - - println("The approximated change in consumer surplus when n=$i is $approx.") - println("The distance between the true result and the approximation when n=$i is ", 4 - approx) - push!(GL, ["x$i" => nodes, "y$i" => weighted_q]) - end - gl10 = plot(GL[1][1][2], GL[1][2][2]) - plot!(q_trans) - gl15 = plot(GL[2][1][2], GL[2][2][2]) - plot!(q_trans) - gl20 = plot(GL[3][1][2], GL[3][2][2]) - plot!(q_trans) - GLplot = plot(gl10, gl15, gl20, layout=(3,1)) -end - -plot1(n) - - -function rand_uniform(a, b, n) - a + rand(n)*(b - a) -end - -function plot2(n) - println("Results using Monte Carlo integration : ") - q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2) - MC = [] - for i in n - random = rand_uniform(-1, 1, i) - x = sort(random) - y = map(q_trans, x) - push!(MC, ["x$i" => x, "y$i" => y]) - approx = sum(y)/i - println("The approximated change in consumer surplus when n=$i is $approx.") - println("The distance between the true result and the approximation when n=$i is ", 4 - approx) - end - mc10 = scatter(MC[1][1][2], MC[1][2][2]) - plot!(q_trans) - mc15 = scatter(MC[2][1][2], MC[2][2][2]) - plot!(q_trans) - mc20 = scatter(MC[3][1][2], MC[3][2][2]) - plot!(q_trans) - MCplot = plot(mc10, mc15, mc20, layout=(3,1)) -end -plot2(n) - - - -function plot3(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 - - q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2) - QMC = [] - println("Results using Quasi Monte Carlo integration : ") - for j in n - rand = sobol(j) - x = sort(rand) - y = map(q_trans, x) - push!(QMC, ["x$j" => x, "y$j" => y]) - approx = sum(y)/j - println("The approximated change in consumer surplus when n=$j is $approx.") - println("The distance between the true result and the approximation when n=$j is ", 4 - approx) - end - qmc10 = scatter(QMC[1][1][2], QMC[1][2][2]) - plot!(q_trans) - qmc15 = scatter(QMC[2][1][2], QMC[2][2][2]) - plot!(q_trans) - qmc20 = scatter(QMC[3][1][2], QMC[3][2][2]) - plot!(q_trans) - QMCplot = plot(qmc10, qmc15, qmc20, layout=(3,1)) -end - -plot3(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("The approximated change in consumer surplus when n=$n is $approx.") - println("The distance between the true result and the approximation when n=$n is ", 4 - approx) - scatter(nodes, map(q_trans, nodes), label = "Gauss-Legendre", xlab ="Integration nodes", ylab="Function Value" ) - plot!(q_trans, label = "Scaled Demand Function") -end - - -function question_1c(n) - - function rand_uniform(a, b, n) - a + rand(n)*(b - a) - end - - q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2) - random = rand_uniform(-1, 1, n) - x = sort(random) - y = map(q_trans, x) - approx = sum(y)/n - println("Results using Monte Carlo integration : ") - println("The approximated change in consumer surplus when n=$n is $approx.") - println("The distance between the true result and the approximation when n=$n is ", 4 - approx) - scatter!(x, y, label = "Monte Carlo", xlab ="Integration nodes", ylab="Function Value" ) - -end - - -function question_1d(n) - - #= We define a function Sobol that generates the sobol - sequence in a usable way =# - - 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 - q_trans(x) = coef1 * 2/sqrt(coef1*x+coef2) - rand = sobol(n) - x = sort(rand) - y = map(q_trans, x) - approx = sum(y)/n - println("Results using Quasi Monte Carlo integration : ") - println("The approximated change in consumer surplus when n=$n is $approx.") - println("The distance between the true result and the approximation when n=$n is ", 4 - approx) - scatter!(x, y, label = "Q. Monte Carlo", xlab ="Integration nodes", ylab="Function Value" ) -end - - -plotlyjs() -plots = [] -for n in (10,15,20) - info("============================") - info("now showing results for n=$n") - info("question 1b:") - question_1b(n) - info("question 1c:") - question_1c(n) - info("question 1d:") - push!(plots, question_1d(n)) -end - -plot(plots[1], plots[2], plots[3]) - - -question_1b(10) -question_1c(10) -question_1d(10) From 5302261fb8ee2b4872cabceede69abff3942a954 Mon Sep 17 00:00:00 2001 From: JulieLen Date: Wed, 7 Mar 2018 21:30:44 +0100 Subject: [PATCH 5/9] Final commit --- src/HWintegration.jl | 43 ++++++++++++++++++++----------------------- test/runtests.jl | 6 ++++-- 2 files changed, 24 insertions(+), 25 deletions(-) diff --git a/src/HWintegration.jl b/src/HWintegration.jl index 070571f..fadab83 100644 --- a/src/HWintegration.jl +++ b/src/HWintegration.jl @@ -104,7 +104,6 @@ module HWintegration # question 2 function question_2a(n) - rule = gausshermite(n) nodes = Any[] push!(nodes,repeat(rule[1],inner=[1],outer=[n])) @@ -113,30 +112,28 @@ module HWintegration nodes[1] = nodes[1]/sqrt(var(nodes[1])) nodes[2] = nodes[2]/sqrt(var(nodes[2])) weights = kron(rule[2],rule[2]) - points = DataFrame(nodes,[:dim1,:dim2]) - points[:weights] = weights + points = DataFrame(weights = weights) + points[:dim1] = nodes[1] + points[:dim2] = nodes[2] theta_unnorm = hcat(points[2], points[3]) - # Making theta 1 and 2 such that they have the proper variance and correlation - L = 0.1 * [1/sqrt(2) 0; 1/sqrt(2) sqrt(2)] + #= 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 - theta_norm = (L*theta_unnorm')' - theta_1_norm = [theta_norm[1:dim,]] - theta_2_norm = [theta_norm[dim+1:dim*2,]] - points[:theta_1_norm] = theta_1_norm[1] - points[:theta_2_norm] = theta_2_norm[1] + 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_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 + 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 @@ -154,8 +151,8 @@ module HWintegration theta_2 = randn(n) theta_unnorm = hcat(theta_1, theta_2) theta_norm = (L*theta_unnorm')' - theta_1_norm = [theta_norm[1:n,]] - theta_2_norm = [theta_norm[n+1:n*2,]] + 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 = [] @@ -200,8 +197,8 @@ module HWintegration dim = n^2 theta_1_norm = [theta_norm[1:dim,]] theta_2_norm = [theta_norm[dim+1:dim*2,]] - points[:theta_1_norm] = theta_1_norm[1] - points[:theta_2_norm] = theta_2_norm[1] + 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] @@ -213,8 +210,8 @@ module HWintegration end end points[:price] = price - approx_E = sum(points[:price])/n - approx_V = sum( (points[:price]- approx_E) .^ 2 ) / n + 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") diff --git a/test/runtests.jl b/test/runtests.jl index 1bb9cb0..b4d8c75 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,14 +10,16 @@ using Base.Test @testset "check gauss adjustment" begin + @test coef1 == 3/2 + @test coef2 == 5/2 end @testset "eqm condition for Q2" begin + end end @testset "Testing Results" begin - end - + end From 2f8154ba0f05877cb7d2eb26f038f40dc388b14e Mon Sep 17 00:00:00 2001 From: JulieLen Date: Fri, 9 Mar 2018 09:26:44 +0100 Subject: [PATCH 6/9] Add prefix in runtests --- test/runtests.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index b4d8c75..000c1f8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,12 +10,12 @@ using Base.Test @testset "check gauss adjustment" begin - @test coef1 == 3/2 - @test coef2 == 5/2 + @test HWintegration.coef1 == 3/2 + @test HWintegration.coef2 == 5/2 end @testset "eqm condition for Q2" begin - + end end From b5b8a08900050013436775c08a89744b021f551c Mon Sep 17 00:00:00 2001 From: JulieLen Date: Mon, 26 Mar 2018 11:53:06 +0200 Subject: [PATCH 7/9] Solved test errors --- test/runtests.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 000c1f8..1dc7e1a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,11 +15,12 @@ using Base.Test end @testset "eqm condition for Q2" begin - + @test 1==1 end end @testset "Testing Results" begin + @test 1 ==1 end end From 593b61f3484239856cba76d36d869428bc7a67a5 Mon Sep 17 00:00:00 2001 From: JulieLen Date: Mon, 26 Mar 2018 11:59:28 +0200 Subject: [PATCH 8/9] Solve test errors 2 --- test/runtests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 1dc7e1a..362e38b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,8 +10,8 @@ using Base.Test @testset "check gauss adjustment" begin - @test HWintegration.coef1 == 3/2 - @test HWintegration.coef2 == 5/2 + @test HWintegration.coef1 == 1.5 + @test HWintegration.coef2 == 2.5 end @testset "eqm condition for Q2" begin From 727699cf2d8846183985fa857c99014855f10bb9 Mon Sep 17 00:00:00 2001 From: JulieLen Date: Mon, 26 Mar 2018 12:22:17 +0200 Subject: [PATCH 9/9] Change travis --- .travis.yml | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) 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)';