# MoneyCycles.R
#
# Copyright (C) 2009 Andrew Clausen <clausen@sas.upenn.edu> and Carlo Strub
# <cs@carlostrub.ch>
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.  It is available here: <http://www.gnu.org/licenses/>.

# load sets library which allows for testing all possible sequences
library(sets)

# choose tolerance for finding roots
tol <- 1e-6

# define exogenous variables
alpha <- 0.3
beta <- 0.98
gamma <- 1.5
cost.bar <- .1
pi <- 0.01

# define functions
u <- function(x) x^alpha
u.inverse <- function(u) u^(1/alpha)
u_ <- function(x) ifelse(x < 0, Inf, alpha * x^(alpha-1))
u_.inverse <- function(u) (u/alpha)^(1/(alpha - 1))

cost.upper <- function(q) cost.bar + q^gamma
cost <- function(q) ifelse(q == 0, 0, cost.upper(q))
cost_ <- function(q) gamma * q^(gamma - 1)
cost_.inverse <- function(cost) (cost/gamma)^(1/(gamma-1))

# net-present marginal utility per period (from the Euler equation)
u_.equilibrium <- function(startx,t) (beta/(1+pi))^(-t) * u_(startx)

# consumption per period
x.equilibrium <- function(startx,t) u_.inverse(u_.equilibrium(startx,t))

# (shadow) marginal cost from first-order conditions 
cost_.shadow <- function(startx,t) u_.equilibrium(startx,t)

# (shadow) production quantity
q.shadow <- function(startx,t) cost_.inverse(cost_.shadow(startx,t))

# inverse of q.shadow for t=0. (given startq, compute startx)
q.shadow.inv <- function(startq) u_.inverse(cost_(startq))

# amount of money left over at the end of the period
budget.m <- function(m, x, q, T) (m + q + T - x) / (1 + pi)

# creates vector of money holdings
budget.m.given.plan <- function(x, q, T)
{
	n <- length(x)
	m <- numeric(n+1)
	m[[1]] <- 0
	for (t in 1:n)
		m[[t+1]] <- budget.m(m[[t]], x[[t]], q[[t]], T)
	m[2:(n+1)]
}

# calculate bounds for maximum plan length.  (See the appendix for an
# explanation.)
calc.bounds <- function(T)
{
	min.startx <- T

	# from proof of Theorem 4
	psi <- function(x) u_(x) - cost_(x - T)
	max.startx <- uniroot(psi, c(min.startx, 100), tol=tol)$root

	# since u' = c', can relate bounds on x to bounds on q.
	min.startq <- q.shadow(max.startx, 0)
	max.startq <- q.shadow(min.startx, 0)

	# from Appendix A
	max.n <- ceiling(log(u_(max.startx) / u_(T)) / log(beta / (1+pi)))

	list(min.startx=min.startx, max.startx=max.startx,
	     min.startq=min.startq, max.startq=max.startq,
	     max.n=max.n, T=T)
}

# returns true if work plan satisfies all budget constraints
plan.is.feasible <- function(bounds, startx) function(workplan)
{
	n <- length(workplan)
	times <- 0:(n-1)
	x <- x.equilibrium(startx, times)
	q <- workplan * q.shadow(startx, times)
	all(budget.m.given.plan(x, q, T) >= 0)
}

# returns value of greedy plan
greedy.plan.value <- function(bounds, n) function(startx)
{
	times <- 0:(n-1)
	plan <- numeric(n)

	V <- 0
	m <- 0
	for (t in times)
	{
		m <- m / (1 + pi) + bounds$T
		x <- x.equilibrium(startx, t)
		if (m < x)
		{
			q <- q.shadow(startx, t)
			V <- V - beta^t * cost(q)
			m <- m + q
			plan[t+1] <- 1
		}
		if (m < x)
			return(-Inf)
		V <- V + beta^t * u(x)
		m <- m - x
	}

	result <- V / (1 - beta^n)
	attr(result, "plan") <- plan
	result
}

# finds optimal startx and its value for a plan of length n
optimal.x.given.n <- function(bounds, n)
{
	result <- optimize(f = greedy.plan.value(bounds, n),
			   lower = bounds$min.startx,
			   upper = bounds$max.startx,
			   maximum = TRUE, tol=tol)
	startx <- result$maximum
	plan <- attr(greedy.plan.value(bounds, n)(startx), "plan")

	times <- 0:(n-1)
	x <- x.equilibrium(startx, times)
	q <- plan * q.shadow(startx, times)
	m_ <- budget.m.given.plan(x, q, bounds$T)
	m <- c(0, m_[1:(n-1)])
	list(x=x, q=q, n=n, m=m, m_=m_, value=result$objective)
}

# finds plan and startx with highest value.
optimal.x <- function(T)
{
	bounds <- calc.bounds(T)

	best <- list(value=-Inf)
	for (n in 1:bounds$max.n)
	{
		this <- optimal.x.given.n(bounds, n)
		if (this$value > best$value)
			best <- this
	}
	best
}

# excess demand given best plan
excess.demand <- function(T)
{
	best <- optimal.x(T)
	(sum(best$x) - sum(best$q)) / best$n
}

# calculate equilibrium by adjusting prices such that markets clear
calc.eq <- function()
{
	T <- uniroot(excess.demand, c(0.00001, 0.1), tol=tol)$root
	eq <- optimal.x(T)
	eq$T <- T

	# measure 1/n agents at each stage of the plan; should equal eq$T / pi
	eq$M <- sum(eq$m_) / eq$n

	eq
}

# write the data out ready for LaTeX to use
write.latex.data <- function(eq)
{
	# create psTricks vectors
	transform.t <- function(t) t * 90/(2*eq$n-1) + 5
	transform.x <- function(x) x * 90/range(eq$q)[2] + 5
	transform.m <- function(x) x * 85/range(eq$m_)[2] + 5
	format.pstricks <- function(x, y, name)
	{
		points <- cbind(x, y)
		points.pairs <- apply(
			points,
			1,
			function(x) paste("{", x[1], ",", x[2], "}"))

		paste("\\savedata{\\", name, "}[{",
		      paste(points.pairs, collapse=","), "}]\n",
		      sep="")
	}

	times <- 1:(2*eq$n) -1
	format.timeseries <- function(x, name)
		format.pstricks(transform.t(times), transform.x(x), name)

	cat(format.timeseries(eq$x, "cyclex"), file="graph")
	cat(format.timeseries(eq$q, "cycleq"), file="graph", append=TRUE)
	cat(format.timeseries(eq$m_, "cyclem"), file="graph", append=TRUE)

	cat(format.pstricks(transform.m(sort(eq$m_)), rep(60, eq$n),
			    "moneyhist"), file="graph", append=TRUE)
}

eq <- calc.eq()
write.latex.data(eq)

# plot best plan
plot(rep(eq$x, 2), type="b", col="blue", ylim=c(0, max(eq$q)),
     xlab="t", ylab="")
lines(rep(eq$m_, 2), type="b", col="green")
lines(rep(eq$q, 2), type="b", col="red")
legend("topright", c("x", "m", "q"), lty=c(1, 1, 1),
       col=c("blue", "green", "red"))