Pre rewrite

ZacLN · Oct 18, 2015 · 34e9b98 · 34e9b98
1 parent a1e462a
commit 34e9b98
Show file tree

Hide file tree

Showing 6 changed files with 96 additions and 59 deletions.
diff --git a/src/EconModel.jl b/src/EconModel.jl
@@ -49,6 +49,7 @@ export  ARSim,
         ∫,
         load,
         updateaggregate!,
+        updateaggregatevariables!,
         updatedistribution!,
         updatetransition!
 

diff --git a/src/aggregatefuncs.jl b/src/aggregatefuncs.jl
@@ -32,7 +32,7 @@ function updatetransition!(M::Model)
         brackets = ntuple(ip->findbracket(Pf[ip][i],M.aggregate.g[ip]),M.state.nendo)
         ei = [findfirst(M.aggregate.G[M.state.nendo+ie][i].==M.aggregate.g[M.state.nendo+ie]) for ie = 1:M.state.nexog]
         abrackets = kron(brackets...)
-        w = spzeros(ntuple(i->length(M.aggregate.g[i]),M.state.nendo)...)
+        w = zeros(ntuple(i->length(M.aggregate.g[i]),M.state.nendo)...)
         for b in abrackets
             w[Tuple(b.i)...]+=b.w
         end
@@ -101,7 +101,7 @@ end
 # end
 
 
-function updateaggregate!(M::Model,ϕ=0.0)
+function updateaggregatevariables!(M::Model,ϕ=0.0)
     if any(M.aggregate.isag)
         for i = 1:M.aggregate.n
             v = M.aggregate.target[i]
@@ -131,3 +131,10 @@ function setaggregate!(M::Model,V::Symbol,ag::Vector{Float64})
         M[V,0,id]=ag[ii]*ones(sum(id))
     end
 end
+
+function updateaggregate!(M,Φ=0.0)
+    updatetransition!(M)
+    updatedistribution!(M)
+    updateaggregatevariables!(M,Φ)
+    return nothing
+end
diff --git a/src/debug.jl b/src/debug.jl
@@ -4,34 +4,62 @@ gtype=CurtisClenshaw
 
 (foc,endogenous,exogenous,policy,static,params,aux,agg)=
 (:[
-  1.0   -  β*R*Expect(Uc[+1]/π[+1])/Uc
-  (1-ϵ)+ϵ*MC-ϕ*π*(π-Sπ)+ϕ/R*Expect(Y[+1]*π[+1]^2*(π[+1]-Sπ))/Y
-],:[],:[
-  ηz      = (1.,0.9,0.001,1)
-  ηr      = (1.,0.9,0.001,0)
+            Π+W*h*η+R/π*b[-1]-b-c
+            (λ*W*η+Uh)*h
+		    (b-blb)-R*β*Expect(1/π[+1]*λ[+1]*(b-blb))/(λ)
+            Γ1  -  (ε*rp^(-ε)*Y*MC  +  θ/R*π^-ε*rp^(-ε)*Expect((1.0/π[+1])^-ε*(1.0/rp[+1])^(-ε)*Γ1[+1]*π[+1]))
+            Γ2  -  ((ε-1.0)*rp^-ε*Y     +  θ/R*π^(1.0-ε)*rp^-ε*Expect((1.0/rp[+1])^-ε*(1.0/π[+1])^(1.0-ε)*Γ2[+1]*π[+1]))
+            Δ   -  (1.0-θ)*rp^-ε     -  θ*(π/π[-1])^ε*Δ[-1]
+            π   - π[-1]/((  (1.0   -  (1.0-θ)*rp^(1.0-ε))  /  θ  )^(1.0/(1.0-ε)))
+
+],:[
+    b       = (-2.,10.,6)
+    Δ       = (1.0,1.001,1)
+    π       = (0.99,1.01,1)
+],:[
+    η       = (1,0.9,0.1,1)
+    ηz       = (1,0.9,0.001,1)
+],:[
+    b       = (-2.,25.,b,0.9)
+    c       = (0,5,0.4)
+    h       = (0,1,.95)
+    π       = (0.99,1.01,1.0)
+    Δ       = (1.0,1.001,1.0)
+    Γ1      = (1,50,12)
+    Γ2      = (1,50,12)
+],:[
+    R       = SR*(π/Sπ)^ϕπ
+    λ 	    = c^-σc
+    Uh  	= -ϕh*(1-h)^-σh
+    hh      = η*h
+    bp      = b*1
+    db      = b-R*b[-1]/π
+    rp      = Γ1/Γ2
+    MC      = W/ηz
+    Y       = ηz*H/Δ
+    Π       = Y-W*H
+    AG      = Y-C-dB
 ],:[
-  C       = (0.1,.5,0.33)
-  π       = (0.97,1.03,1)
+    β       = 0.98
+    σc      = 2.5
+    ϕh      = 6.5
+    σh      = 2.0
+    blb     = -2.0
+    ε       = 10.0
+    θ       = 0.75
+    ϕπ      = 2.5
+    ϕY      = 0.1
+    Sπ      = 1.00
+    SR      = 1.015
 ],:[
-  Y       = C/(1.0-ϕ/2.0*(π-Sπ)^2.0)
-  H       =  Y/ηz
-  Uc      = C^-σc
-  Uh      = -ϕh*(1.0-H)^-σh
-  W       = -Uh/Uc
-  MC      = W/ηz
-  Z       = SR*((π/Sπ)^ϕπ)*ηr
-  R       = (-(((-1.0)+(-Z))-(((-1.0)-(-Z))^2.0)^.5)/2.0)
+    W = (0.9,1.0)
 ],:[
-  β       = 0.98
-  σc      = 2.5
-  σh      = 2.0
-  ϕh      = 6.5
-  ϵ       = 10.0
-  ϕ       = 110.0
-  ϕπ      = 2.5
-  Sπ      = 1.00
-  SR      = $(1.00/0.98)
-],:[],:[])
+    ηz
+    B   = (b,0)
+    C   = (c,0.3)
+    H   = (hh,0.3)
+    dB  = (db,0)
+])
 
 
 

diff --git a/src/model.jl b/src/model.jl
@@ -1,27 +1,26 @@
 type ModelMeta
-  foc::Expr
-  parameters::Expr
-  endogenous::Expr
-  exogenous::Expr
-  policy::Expr
-  static::Expr
-  auxillary::Expr
-  aggregate::Expr
+    foc::Expr
+    parameters::Expr
+    endogenous::Expr
+    exogenous::Expr
+    policy::Expr
+    static::Expr
+    auxillary::Expr
+    aggregate::Expr
 end
 
 type Model
-  F::Function
-  J::Function
-  E::Function
-  state::StateVariables
-  policy::PolicyVariables
-  future::FutureVariables
-  static::StaticVariables
-  auxillary::AuxillaryVariables
-  aggregate::AggregateVariables
-  error::Array{Float64,2}
-  meta::ModelMeta
-  # D
+    F::Function
+    J::Function
+    E::Function
+    state::StateVariables
+    policy::PolicyVariables
+    future::FutureVariables
+    static::StaticVariables
+    auxillary::AuxillaryVariables
+    aggregate::AggregateVariables
+    error::Array{Float64,2}
+    meta::ModelMeta
 end
 
 
@@ -40,7 +39,7 @@ function Model(foc::Expr,endogenous::Expr,exogenous::Expr,policy::Expr,static::E
     @assert length(foc.args) == length(policy.args) "equations doesn't equal numer of policy variables"
 
 
-    meta                   = ModelMeta(deepcopy(foc),
+    meta                    = ModelMeta(deepcopy(foc),
                                         deepcopy(params),
                                         deepcopy(endogenous),
                                         deepcopy(exogenous),
@@ -86,12 +85,11 @@ function Model(foc::Expr,endogenous::Expr,exogenous::Expr,policy::Expr,static::E
         if !in(aux.args[i].args[1],[x.args[1] for x in variablelist[:,1]])
             warn("Added $(aux.args[i].args[1]) to variable list") # this may no longer be hit ever
             x = copy(aux.args[i].args[1])
-            # x = addindex!(x,iglist)
             x = addindex!(x)
             x = hcat(x,:(M.auxillary.X[i,$i]),symbol("A$i"))
             variablelist = vcat(variablelist,x)
         end
-    end 
+    end
 
     for i = State.nendo+1:State.n
         if !in(State.names[i],[x.args[1] for x in variablelist[:,1]])

diff --git a/src/optim.jl b/src/optim.jl
@@ -116,20 +116,24 @@ function ridders(f::Function, a::Real, b::Real;
 	while niter < maxiter && abs(x1 - x2) > tol
 		xm = (x1 + x2)/2.0
 		fm = f(xm)
-		if fm == 0; return xm; end
+		abs(fm) < tol ? (return xm) : nothing
 
 		x3 = xm + (xm - x1) * sign(f1 - f2) * fm / sqrt(fm^2 - f1 * f2)
 		f3 = f(x3)
 		niter += 2
-		if f3 == 0; return x3; end
+		abs(f3) < tol ? (return x3) : nothing
 
 		if fm * f3 < 0
-			x1 = xm;  f1 = fm
-			x2 = x3;  f2 = f3
+			x1 = xm
+            f1 = fm
+			x2 = x3
+            f2 = f3
 		elseif f1 * f3 < 0
-			x2 = x3;  f2 = f3
+			x2 = x3
+            f2 = f3
 		elseif f2 * f3 < 0
-			x1 = x3;  f1 = f3
+			x1 = x3
+            f1 = f3
 		else
 			error("Inform the maintainer: you should never get here.")
 		end

diff --git a/src/solve.jl b/src/solve.jl
@@ -59,7 +59,7 @@ function solve!(M::Model,
 
     for iter = 1:n
         maximum(abs(M.error))<crit*10 ? upf = 1 : nothing
-        (mod(iter,upag)==0 || maximum(abs(M.error))<crit) && M.aggregate.n>0  &&  upag != -1 ? (updatetransition!(M);updatedistribution!(M);updateaggregate!(M,Φ)) : nothing
+        (mod(iter,upag)==0 || maximum(abs(M.error))<crit) && M.aggregate.n>0  &&  upag != -1 ? (updatetransition!(M);updatedistribution!(M);updateaggregatevariables!(M,Φ)) : nothing
         (mod(iter,upf)==0 ||iter <4) ? (getfuture(M);updatefutureauxillary(M);updatefutureaggregate(M)) : nothing
 
         M.E(M)
@@ -72,8 +72,7 @@ function solve!(M::Model,
 
         for i = 1:M.state.G.n
             x = vec(M.policy.X[i,:])-vec(M.J(M,i)\vec(M.error[i,:]))
-            x = max(M.policy.lb,x)
-            x = min(M.policy.ub,x)
+            x = min(M.policy.ub,max(M.policy.lb,x))
             @inbounds M.policy.X[i,:] = vec(M.policy.X[i,:])*ϕ +  (1-ϕ)*x
         end