Commit Daniel and Luis HW

src/HWintegration_LB_DB.jl

@@ -0,0 +1,146 @@
+module HWintegration
+
+	const A_SOL = 4
+	println("Analytically, the change of CR equals 4")
+	# imports
+	using FastGaussQuadrature
+	using Roots
+	using Sobol
+	using Plots
+	using Distributions
+	using NullableArrays
+	using DataFrames
+	using DataStructures  # OrderedDict
+
+	# set random seed
+	srand(12345)
+
+	# demand function
+	function f(x)
+		2*x^(-.5)
+	end
+
+	# makes a plot for question 1
+	function plot_q1()
+	global dem = plot(f, 0.1, 10, title="Demand", xaxis=("p"), yaxis=("quantitiy"), labels = "f(x)")
+	hline!([2, 1], label = "")
+	end
+
+	function question_1a(np)
+		for n in (10,15,20)
+			np = gausslegendre(n)
+			glx = 3/2*np[1] + 3/2
+			gly = [f(3/2*x + 5/2) for x in np[1]]
+			integ = 3/2*sum([f(3/2*x + 5/2) for x in np[1]] .* np[2])
+			global gl = scatter(glx, gly, labels = "GL", title = "$n nodes")
+			error = A_SOL-integ
+			println("GL-M: integral = $integ and distance to optimal sol. = $error (with $n nodes)")
+	end
+end
+
+
+	function question_1b(n)
+		for n in (10,15,20)
+			mcx = rand(Uniform(1,4), n)
+			mcy = f.(mcx)
+			integ = 3/n * sum([f(x) for x in mcx])
+			global mc = scatter(mcx, mcy, labels = "MC", title = "$n nodes")
+			error = A_SOL-integ
+			println("MC: integral = $integ and distance to optimal sol. = $error (with $n nodes)")
+	end
+end
+
+	function question_1c(n)
+		for n in (10,15,20)
+			s = SobolSeq(1)
+			qmcx = 3*[ hcat([next(s) for i = 1:n])[x][1] for x = 1:n]+1
+			qmcy = f.(qmcx)
+			integ = 3/n * sum([f(x) for x in qmcx])
+			global qmc = scatter(qmcx, qmcy, labels = "Quasi MC", title = "$n nodes")
+			println("Quasi MC: integral = $integ and distance to optimal sol. = $error (with $n nodes)")
+	end
+end
+
+	# question 2
+# theta1*p^(-1) + theta2*p^(-0.5) = 2
+
+	function question_2a()
+		function p(t = [0,0])
+			function eq(x)
+				exp(t[1])* float(x)^(-1) + exp(t[2])*float(x)^(-0.5) - 2
+			end
+			fzero(eq, [0,1000])
+		end
+		immutable MvNormal{Cov<:AbstractPDMat,Mean<:Vector} <: AbstractMvNormal
+		            μ::Mean
+		            Σ::Cov
+		end
+		μ = [0,0]
+		Σ = [0.02   0.01 ; 0.01  0.01]
+		rand(MvNormal(μ, Σ), n)
+
+		nodes = Any[]
+		push!(nodes,repeat(rules["hermite"][1],inner=[1],outer=[10]))  # dim1
+		push!(nodes,repeat(rules["hermite"][1],inner=[2],outer=[5]))  # dim2
+		weights = kron(rules["hermite"][2],kron(rules["hermite"][2],rules["hermite"][2]))
+		df = hcat(DataFrame(weights=weights),DataFrame(nodes,[:dim1,:dim2]))
+
+	end
+
+	function question_2b(n)
+		for n in (10,15,20)
+		μ = [0,0]
+		Σ = [0.02   0.01 ; 0.01  0.01]
+		draws = rand(MvNormal(μ, Σ),  n)
+		draws = [[draws[1,i], draws[2,i]] for i in 1:n]
+		function p(t = [0,0])
+			function eq(x)
+				exp(t[1])* float(x)^(-1) + exp(t[2])*float(x)^(-0.5) - 2
+			end
+			fzero(eq, [0,1000])
+		end
+		expec= 1/n * sum(p.(draws))
+
+		var = 1/n * sum((p.(draws) - expec).^2)
+		println("By Monte Carlo, the expectation is $expec and the variance is $var")
+	end
+end
+
+	function question_2bonus()
+
+	end
+
+	# function to run all questions
+	function runall()
+		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_1a(n)
+			info("question 1c:")
+			question_1b(n)
+			info("question 1d:")
+			question_1c(n)
+			println("")
+			allplots = plot(dem, gl, mc, qmc, layout = (2,2))
+			display(allplots)
+			savefig("n=$n.png")
+		end
+			info("question 2a:")
+			q2 = question_2a(n)
+			println(q2)
+			info("question 2b:")
+			for n in (10,15,20)
+			info("Now showing results for n = $n")
+			question_2b(n)
+			println(q2b)
+		end
+			info("bonus question: Quasi monte carlo:")
+			q2bo = question_2bonus(n)
+			println(q2bo)
+		end
+	info("end of HWintegration")
+end