it was very frustrating to find out that the java portable function that I wrote for rate() function from VBA financial was wrong.
Finally after 2 days of painstaking search and lot of hours in front of monitor.. I could get the code of java based rate() function together.
I don't understand ABCD of this function's internals. Only thing that I know is it works :-) and I guess thats what matters too.
I got lot of help from the GNUcash people as well as from Gnumeric and openoffice implementations.
I hope this function works rightly.
Hope it will be useful for someone.
Use at your own liability but can do anything with the source code otherwise.
/*
* Created on Mar 10, 2005
*/
package test;
/**
* @author vijay_dharap This is rate derived from source code of Gnumeric
*
*/
public class GnumericRate {
private static int GOAL_SEEK_OK = 0;
private static int GOAL_SEEK_ERROR = 1;
public static void main(String[] args) {
GnumericRate rateCalculator = new GnumericRate();
// double rate = gnumeric_rate(36,-38035, 4910, 21252,1);
// double rate = gnumeric_rate(36,-38035, 1400, 21252,1);
// double rate = gnumeric_rate(36,-38035, 3781.2, 21252,1);
// double rate = gnumeric_rate(36, -38035, 8500, 21252, 1);
double rate = rateCalculator.rate(36, -38035, 18500, 21252, 1);
System.out.println(rate);
}
private double gnumeric_rate_f(double rate/* 0.1 */, double nper,
double pv, double pmt, double fv, int type) {
return pv * calculate_pvif(rate, nper) + pmt * (1 + rate * type)
* calculate_fvifa(rate, nper) + fv;
}
/* The derivative of the above function with respect to rate. */
private double gnumeric_rate_df(double rate/* 0.1 */, double nper,
double pv, double pmt, int type) {
return -pmt * calculate_fvifa(rate, nper) / rate
+ calculate_pvif(rate, nper - 1) * nper
* (pv + pmt * (type + 1 / rate));
}
private double rate(/* double rate, */double nper, double pv,
double pmt, double fv, int type) {
GoalSeekData data = new GoalSeekData();
double rate0 = 0.1; // hard coding
int status; // used for deciding if NR method was successful
data.xmin = Math.max(data.xmin, -Math.pow(Double.MAX_VALUE / 1e10,
1.0 / nper) + 1);
data.xmax = Math.min(data.xmax, Math.pow(Double.MAX_VALUE / 1e10,
1.0 / nper) - 1);
/* Newton search from guess. */
status = goal_seek_newton(data, rate0, nper, pv, pmt, fv, type);
System.out.println("0--> ok, 1-->error: " + status);
if (status != GOAL_SEEK_OK) {
int factor;
/* Lay a net of test points around the guess. */
for (factor = 2; !(data.havexneg && data.havexpos) && factor < 100; factor *= 2) {
goal_seek_point( /* data, */rate0 * factor, nper, pv, pmt, fv,
type);
goal_seek_point( /* data, */rate0 / factor, nper, pv, pmt, fv,
type);
}
/* Pray we got both sides of the root. */
status = goal_seek_bisection(data, nper, pv, pmt, fv, type);
System.out.println("0--> ok, 1-->error: " + status);
}
return data.root;
}
/*
* Seek a goal using a single point.
*/
private double goal_seek_point(/* GoalSeekData data, */
double x0/* 0.1 */, double nper, double pv, double pmt, double fv, int type) {
double y0;
y0 = gnumeric_rate_f(x0, /* &y0, */nper, pv, pmt, fv, type);
// as far my understanding of code goes, this is not required as
// we are not botheredwith what update_datareturns
// y0 = update_data (x0, y0, data);
return y0;
}
/*
* Seek a goal (root) using Newton's iterative method.
*
* The supplied function must (should) be continously differentiable in the
* supplied interval. If NULL is used for `df', this function will estimate
* the derivative.
*
* This method will find a root rapidly provided the initial guess, x0, is
* sufficiently close to the root. (The number of significant digits
* (asympotically) goes like i^2 unless the root is a multiple root in which
* case it is only like c*i.)
*/
private int goal_seek_newton(GoalSeekData data, double x0/* 0.1 */,
double nper, double pv, double pmt, double fv, int type) {
double precision = data.precision / 2;
// as of now not givnig null. to test how it works
Object df = new Object();
System.out.println("goal_seek_newton");
for (int iterations = 0; iterations < 20; iterations++) {
double x1, y0, df0, stepsize;
// status is to be removed
y0 = gnumeric_rate_f(x0, /* &y0, */nper, pv, pmt, fv, type);
System.out.println("x0 = " + x0);
System.out.println("y0 = " + y0);
if (update_data(x0, y0, data))
return GOAL_SEEK_OK;
if (df != null)
df0 = gnumeric_rate_df(x0, /* &df0, */nper, pv, pmt, type);
else {
double xstep;
if (Math.abs(x0) < 1e-10)
if (data.havexneg && data.havexpos)
xstep = Math.abs(data.xpos - data.xneg) / 1e6;
else
xstep = (data.xmax - data.xmin) / 1e6;
else
xstep = Math.abs(x0) / 1e6;
df0 = fake_df(x0, /* &df0, */xstep, data, nper, pv, pmt, fv,
type);
}
/* If we hit a flat spot, we are in trouble. */
if (df0 == 0) {
return GOAL_SEEK_ERROR;
}
/*
* Overshoot slightly to prevent us from staying on just one side of
* the root.
*/
x1 = x0 - 1.000001 * y0 / df0;
if (x1 == x0) {
data.root = x0;
return GOAL_SEEK_OK;
}
stepsize = Math.abs(x1 - x0) / (Math.abs(x0) + Math.abs(x1));
System.out.println("df0 = " + df0);
System.out.println("ss = " + stepsize);
x0 = x1;
if (stepsize < precision) {
data.root = x0;
return GOAL_SEEK_OK;
}
}
return GOAL_SEEK_ERROR;
}
/*
* Seek a goal (root) using bisection methods.
*
* The supplied function must (should) be continous over the interval.
*
* Caller must have located a positive and a negative point.
*
* This method will find a root steadily using bisection to narrow the
* interval in which a root lies.
*
* It alternates between mid-point-bisection (semi-slow, but guaranteed
* progress), secant-bisection (usually quite fast, but sometimes gets
* nowhere), and Ridder's Method (usually fast, harder to fool than the
* secant method).
*/
private int goal_seek_bisection(GoalSeekData data, double nper,
double pv, double pmt, double fv, int type) {
final int M_SECANT = 0;
final int M_RIDDER = 1;
final int M_NEWTON = 2;
final int M_MIDPOINT = 3;
int method = 0; // default m_secant
int iterations;
double stepsize;
int newton_submethod = 0;
System.out.println("goal_seek_bisection");
if (!data.havexpos || !data.havexneg)
return GOAL_SEEK_ERROR;
stepsize = Math.abs(data.xpos - data.xneg)
/ (Math.abs(data.xpos) + Math.abs(data.xneg));
/* log_2 (10) = 3.3219 < 4. */
final int GNUM_DIG = 15; // dont know what this is
for (iterations = 0; iterations < 100 + GNUM_DIG * 4; iterations++) {
double xmid, ymid;
int status;
method = (iterations % 4 == 0) ? M_RIDDER
: ((iterations % 4 == 2) ? M_NEWTON : M_MIDPOINT);
//again:
switch (method) {
default:
System.exit(2);
case M_SECANT:
xmid = data.xpos - data.ypos
* ((data.xneg - data.xpos) / (data.yneg - data.ypos));
break;
case M_RIDDER:
double det;
xmid = (data.xpos + data.xneg) / 2;
ymid = gnumeric_rate_f(xmid, /* &ymid, */nper, pv, pmt, fv,
type);
if (ymid == 0) {
update_data(xmid, ymid, data);
return GOAL_SEEK_OK;
}
det = Math.sqrt(ymid * ymid - data.ypos * data.yneg);
if (det == 0)
/* This might happen with underflow, I guess. */
continue;
xmid += (xmid - data.xpos) * ymid / det;
break;
case M_MIDPOINT:
xmid = (data.xpos + data.xneg) / 2;
break;
case M_NEWTON:
double x0,
y0 = 0,
xstep,
df0;
/* This method is only effective close-in. */
if (stepsize > 0.1) {
method = M_MIDPOINT;
xmid = (data.xpos + data.xneg) / 2;
// goto again;
} else {
switch (newton_submethod++ % 4) {
case 0:
x0 = data.xpos;
x0 = data.ypos;
break;
case 2:
x0 = data.xneg;
y0 = data.yneg;
break;
default:
case 3:
case 1:
x0 = (data.xpos + data.xneg) / 2;
y0 = gnumeric_rate_f(x0, /* &y0, */nper, pv, pmt, fv,
type);
}
xstep = Math.abs(data.xpos - data.xneg) / 1e6;
df0 = fake_df(x0, /* &df0, */xstep, data, nper, pv, pmt,
fv, type);
if (df0 == 0)
continue;
/*
* Overshoot by 1% to prevent us from staying on just one
* side of the root.
*/
xmid = x0 - 1.01 * y0 / df0;
if ((xmid < data.xpos && xmid < data.xneg)
|| (xmid > data.xpos && xmid > data.xneg))
/* We left the interval. */
continue;
}
}
ymid = gnumeric_rate_f(xmid, /* &ymid, */nper, pv, pmt, fv, type);
String themethod;
switch (method) {
case M_MIDPOINT:
themethod = "midpoint";
break;
case M_RIDDER:
themethod = "Ridder";
break;
case M_SECANT:
themethod = "secant";
break;
case M_NEWTON:
themethod = "Newton";
break;
default:
themethod = "?";
}
System.out.println("xmid = " + xmid + " themethod= " + themethod);
System.out.println(" ymid = " + ymid);
if (update_data(xmid, ymid, data)) {
return GOAL_SEEK_OK;
}
stepsize = Math.abs(data.xpos - data.xneg)
/ (Math.abs(data.xpos) + Math.abs(data.xneg));
System.out.println("ss = " + stepsize);
if (stepsize < data.precision) {
if (data.yneg < ymid) {
ymid = data.yneg;
xmid = data.xneg;
}
if (data.ypos < ymid) {
ymid = data.ypos;
xmid = data.xpos;
}
data.root = xmid;
return GOAL_SEEK_OK;
}
}
return GOAL_SEEK_ERROR;
}
/* Calculate (1+x)^r accurately. */
private double pow1p(double x, double y) {
if (Math.abs(x) > 0.5)
return Math.pow(1 + x, y);
else
return Math.exp(y * Math.log1p(x));
}
/* Calculate ((1+x)^r)-1 accurately. */
private double pow1pm1(double x, double y) {
if (x <= -1)
return Math.pow(1 + x, y) - 1;
else
return Math.expm1(y * Math.log1p(x));
}
private double calculate_pvif(double rate, double nper) {
return pow1p(rate, nper);
}
private double calculate_fvifa(double rate, double nper) {
/* Removable singularity at rate == 0. */
if (rate == 0)
return nper;
else
return pow1pm1(rate, nper) / rate;
}
/*
* Calculate a reasonable approximation to the derivative of a function in a
* single point.
*/
private double fake_df(double x, /* double *dfx, */double xstep,
GoalSeekData data, double nper, double pv, double pmt, double fv,
int type) {
double dfx;
double xl, xr, yl, yr;
System.out.println("fake_df (x= " + x + ", xstep= " + xstep);
xl = x - xstep;
if (xl < data.xmin)
xl = x;
xr = x + xstep;
if (xr > data.xmax)
xr = x;
if (xl == xr) {
System.out.println("==> xl == xr\n");
return GOAL_SEEK_ERROR;
}
yl = gnumeric_rate_f(xl, /* &yl, */nper, pv, pmt, fv, type);
System.out.println("==> xl= " + xl + " yl= " + yl);
yr = gnumeric_rate_f(xr, /* &yr, */nper, pv, pmt, fv, type);
System.out.println("==> xr= " + xr + " yr= " + yr);
dfx = (yr - yl) / (xr - xl);
System.out.println("==> " + dfx);
return dfx;
}
/*
* I DONT UNDERSTAND ANYTHING AT ALL OF THIS METHOD
*/
private boolean update_data(double x, double y, GoalSeekData data) {
if (y > 0) {
if (data.havexpos) {
if (data.havexneg) {
/*
* When we have pos and neg, prefer the new point only if it
* makes the pos-neg x-internal smaller.
*/
if (Math.abs(x - data.xneg) < Math.abs(data.xpos
- data.xneg)) {
data.xpos = x;
data.ypos = y;
}
} else if (y < data.ypos) {
/* We have pos only and our neg y is closer to zero. */
data.xpos = x;
data.ypos = y;
}
} else {
data.xpos = x;
data.ypos = y;
data.havexpos = true;
}
return false;
} else if (y < 0) {
if (data.havexneg) {
if (data.havexpos) {
/*
* When we have pos and neg, prefer the new point only if it
* makes the pos-neg x-internal smaller.
*/
if (Math.abs(x - data.xpos) < Math.abs(data.xpos
- data.xneg)) {
data.xneg = x;
data.yneg = y;
}
} else if (-y < -data.yneg) {
/* We have neg only and our neg y is closer to zero. */
data.xneg = x;
data.yneg = y;
}
} else {
data.xneg = x;
data.yneg = y;
data.havexneg = true;
}
return false;
} else {
/* Lucky guess... */
data.root = x;
return true;
}
}
}
class GoalSeekData {
double xmin; /* Minimum allowed value for x. */
double xmax; /* Maximum allowed value for x. */
double precision; /* Desired relative precision. */
boolean havexpos; /* Do we have a valid xpos? */
double xpos; /* Value for which f(xpos) > 0. */
double ypos; /* f(xpos). */
boolean havexneg; /* Do we have a valid xneg? */
double xneg; /* Value for which f(xneg) < 0. */
double yneg; /* f(xneg). */
double root; /* Value for which f(root) == 0. */
public GoalSeekData() {
havexpos = havexneg = false;
xmin = -1e10;
xmax = +1e10;
precision = 1e-10;
}
}
any comments are welcome
