Hi everyone,
To solve equation below to calculate the unknown r:
Data looks like:
data test; input stock$ year P B B1 B2 B3 B4 F F1 F2 F3 F4 F5 g; datalines; A 1 80.16 12.1492 14.7642 17.7592 21.1642 24.6292 3.8886 5.23 5.99 6.81 6.93 6.63 0.3038 A 2 63.54 15.3785 18.4785 21.9535 25.8685 29.9185 4.8659 6.2 6.95 7.83 8.1 8.51 0.3769 A 3 41.19 11.5507 15.1507 19.1257 23.5857 28.3907 4.5454 7.2 7.95 8.92 9.61 10.47 0.205 A 4 35.28 15.1777 18.7327 22.7277 27.1727 32.2727 4.3858 7.11 7.99 8.89 10.2 11.14 0.1478 A 5 52 14.9545 18.3545 22.1945 26.4395 31.1195 5.7079 6.8 7.68 8.49 9.36 11.07 0.0644 A 6 50.78 15.8837 19.6137 23.8137 28.4137 33.5537 6.1861 7.46 8.4 9.2 10.28 11.43 -0.0636 A 7 67.64 19.3065 23.2765 27.6065 32.3565 37.6015 6.7226 7.94 8.66 9.5 10.49 12.24 -0.0915 A 8 72.05 21.6415 26.5265 31.9415 37.8015 44.3815 7.7686 9.77 10.83 11.72 13.16 14.59 -0.0553 A 9 74.51 23.6561 30.0461 36.9961 44.5511 52.4911 8.804 12.78 13.9 15.11 15.88 11.32 0.0005 A 10 43.05 27.2373 34.5373 42.7023 51.6873 61.0373 9.0739 14.6 16.33 17.97 18.7 16.621 0.123 B 1 37.8 6.3668 7.6503 9.0298 10.5533 12.2488 2.7893 2.567 2.759 3.047 3.391 3.807 0.3038 B 2 35.4 7.6967 8.9367 10.2797 11.7422 13.3282 3.238 2.48 2.686 2.925 3.172 3.46 0.3769 B 3 29.3 8.6812 10.0817 11.5857 13.2137 14.9497 3.2122 2.801 3.008 3.256 3.472 3.778 0.205 B 4 30.9 8.999 10.3945 11.8675 13.4585 15.18 3.2352 2.791 2.946 3.182 3.443 3.753 0.1478 B 5 29.75 9.5064 10.9409 12.4669 14.1014 15.9644 3.5381 2.869 3.052 3.269 3.726 3.942 0.0644 B 6 39.3 10.1293 11.6418 13.2513 14.9848 16.8748 3.4308 3.025 3.219 3.467 3.78 3.954 -0.0636 B 7 43.3 12.2344 13.8109 15.4949 17.3094 19.3754 3.7561 3.153 3.368 3.629 4.132 4.534 -0.0915 B 8 52 13.2885 14.78 16.4125 18.216 20.2585 4.2523 2.983 3.265 3.607 NA NA -0.0553 B 9 41.6 13.9864 15.6949 17.5539 19.5994 21.7379 4.7464 3.417 3.718 4.091 4.277 4.836 0.0005 B 10 50.2 14.043 16.067 18.2475 20.607 23.178 4.2598 4.048 4.361 4.719 5.142 5.474 0.123 C 1 24.19 1.4639 3.5839 4.8589 6.6539 9.1189 0.3939 4.24 2.55 3.59 4.93 5.63 0.3038 C 2 19.52 2.3384 3.9534 5.7434 7.9884 10.8534 1.3403 3.23 3.58 4.49 5.73 6.51 0.3769 C 3 7.96 2.626 3.711 5.381 7.581 10.846 -1.0251 2.17 3.34 4.4 6.53 5.9 0.205 C 4 15.72 2.7283 4.0083 5.8033 8.0983 11.6433 0.3156 2.56 3.59 4.59 7.09 7.28 0.1478 C 5 11.81 2.3651 3.9451 5.7001 7.9701 10.8651 0.049 3.16 3.51 4.54 5.79 6.39 0.0644 C 6 16.9 2.8771 4.4821 6.5371 9.0671 11.6921 0.964 3.21 4.11 5.06 5.25 7.21 -0.0636 C 7 18.84 3.532 5.792 8.587 12.027 16.417 1.1485 4.52 5.59 6.88 NA NA -0.0915 C 8 21.66 4.5199 7.5199 11.0849 15.4049 20.4099 2.0169 6 7.13 8.64 10.01 11.18 -0.0553 C 9 12.75 4.379 7.449 11.774 16.944 23.244 2.0582 6.14 8.65 10.34 12.6 13.89 0.0005 C 10 24 4.6028 8.4778 13.4228 19.5128 26.4078 0.9412 7.75 9.89 12.18 13.79 16.74 0.123 ;
This is a variant of a previous post with a different data layout (https://communities.sas.com/t5/SAS-Programming/Solving-Equation/m-p/766280#M242843).
Thanks a lot for the sample data. I checked several cases "manually" and found different scenarios:
I've implemented a search for initial values checking the sign of the polynomial function k at r = 0, 0.001, 0.002, ..., 0.999, 1. The midpoint of an interval with a sign change between its boundaries is selected as an initial value (taking into account even the unlikely case that an interval boundary happens to be a root of k). It turned out that, at least for the 483−67=416 usable sample data observations, the initial values found by this algorithm are close enough to the roots of k in (0, 1), if any, so that the SOLVE function finds the roots accurately. The same roots were found with different subdivisions of [0, 1], e.g., into 10000 subintervals (instead of 1000) or only 3 (!). The interval boundaries enclosing a sign change could also be used as starting values for the bisection method. This would be an alternative approach which could replace the SOLVE function if it didn't work in "extreme" cases.
But I've checked all 171 solutions found for the 416 complete observations (163 with a single solution in (0, 1), 4 with two distinct solutions in (0, 1), 249 with no solutions in (0, 1)) with the code used in the earlier post (i.e., using the rational function h rather than the polynomial), except that no renaming of variables was necessary. These checks resulted in a maximum absolute value of h(r) of 6.08E-6. I took a closer look at the case where this relatively "large" maximum occurred (stock='CH0008742519' & year=1998) and found that
The solutions found with the above approach also confirmed those that I had investigated manually.
Here is the code:
/* Define function XRATE2 */
proc fcmp outlib=work.funcs.test;
function k(r, p, B[*], F[*], g);
c6= p;
c5= p*(5-g)-F[1]/2;
c4= p*(10-5*g)-B[1]+B[2]+(g-5)*F[1]/2-F[2]/2;
c3= 10*p*(1-g)+(g-3)*B[1]+(2-g)*B[2]+B[3]+((5*g-9)*F[1]+(g-4)*F[2]-F[3])/2;
c2= p*(5-10*g)+3*(g-1)*B[1]+(1-2*g)*B[2]+(1-g)*B[3]+B[4]+((9*g-7)*F[1]+(4*g-5)*F[2]+(g-3)*F[3]-F[4])/2;
c1= p*(1-5*g)+(3*g-1)*B[1]-g*(B[2]+B[3]+B[4])-(1+g)*B[5]+((7*g-2)*F[1]+(5*g-2)*F[2]+(3*g-2)*F[3]+(g-2)*F[4])/2-F[5];
c0= -p*g+g*B[1]+g*(F[1]+F[2]+F[3]+F[4])-F[5];
z = c6*r**6 + c5*r**5 + c4*r**4 + c3*r**3 + c2*r**2 + c1*r + c0;
return(z);
endfunc;
function xrate2(i, p, B[*], F[*], g); /* first argument is the initial value */
array solvopts[1] i;
r=solve("k", solvopts, 0, ., p, B, F, g);
if min(abs(1+r),abs(r-g))<1e-4 then r=.; /* avoid poles of the original function */
return(r);
endfunc;
run;
/* Make function available */
options cmplib=work.funcs;
/* Find initial values for solutions in (0, 1) */
%let nsubintv=1000; /* number of subintervals of [0, 1] to search */
data init(drop=_:);
set have;
if nmiss(P, g, of F1-F5, of B0-B4) then output; /* Solutions will be missing for incomplete data. */
else do;
array init[6]; /* Up to six roots are possible theoretically. */
array B[*] B0-B4;
array F[*] F1-F5;
do _i=0 to &nsubintv;
_x=_i/&nsubintv;
_z=k(_x, p, B, F, g);
_s=sign(_z);
_prev_s=lag(_s);
if _s~=0 then do;
if _i & _s~=_prev_s~=0 then do;
_j=sum(_j,1);
init[_j]=(_i-0.5)/&nsubintv; /* select midpoint of subinterval with sign change between boundaries */
end;
end;
else do; /* This case (solution on subinterval boundary) is very unlikely. */
_j=sum(_j,1);
init[_j]=_x;
end;
end;
output;
end;
run;
/* Apply function XRATE2 to dataset INIT */
data want(drop=_j);
set init;
array B[*] B0-B4;
array F[*] F1-F5;
array init[6];
array r[6];
do _j=1 to dim(r) while(init[_j]>.);
r[_j]=xrate2(init[_j], p, B, F, g); /* solution found with initial value init[_j] */
end;
run;
The PROC FCMP step and the OPTIONS statement are exactly the same as in my earlier post. I repeated them because your attachment sample.sas contained something different. For your sample data the dimension of the arrays init and r in the DATA steps could be reduced from 6 to 2 (remaining variables are always missing), but theoretically a polynomial of grade 6 can have up to 6 real roots, so it's better to allow for that many solutions and initial values. You may want to drop variables init1-init6. I've left them in dataset WANT just to check the differences (always within ±0.0005 for &nsubintv=1000) between them and the corresponding solutions.
Did you try to modify @FreelanceReinhs solution?
Hello @MAC1430,
@MAC1430 wrote:
yes, I tried and still trying, got error "ERROR 71-185: The xrate function call does not have enough arguments".
Please show your code (in a code box opened with the "running man" button) including both the PROC FCMP step defining the new XRATE function and the DATA step calling it. Then we'll have a chance to diagnose what went wrong. In particular, we'll see how you defined the function and what arguments you may be missing in the function call.
Thank you very much Reinhard, below is my current try (could be silly) by modifying your code.
data test;
input stock$ y P B B1 B2 B3 B4 F F1 F2 F3 F4 F5 g;
datalines;
A 1 80.16 12.1492 14.7642 17.7592 21.1642 24.6292 3.8886 5.23 5.99 6.81 6.93 6.63 0.3038
A 2 63.54 15.3785 18.4785 21.9535 25.8685 29.9185 4.8659 6.2 6.95 7.83 8.1 8.51 0.3769
A 3 41.19 11.5507 15.1507 19.1257 23.5857 28.3907 4.5454 7.2 7.95 8.92 9.61 10.47 0.205
A 4 35.28 15.1777 18.7327 22.7277 27.1727 32.2727 4.3858 7.11 7.99 8.89 10.2 11.14 0.1478
A 5 52 14.9545 18.3545 22.1945 26.4395 31.1195 5.7079 6.8 7.68 8.49 9.36 11.07 0.0644
A 6 50.78 15.8837 19.6137 23.8137 28.4137 33.5537 6.1861 7.46 8.4 9.2 10.28 11.43 -0.0636
A 7 67.64 19.3065 23.2765 27.6065 32.3565 37.6015 6.7226 7.94 8.66 9.5 10.49 12.24 -0.0915
A 8 72.05 21.6415 26.5265 31.9415 37.8015 44.3815 7.7686 9.77 10.83 11.72 13.16 14.59 -0.0553
A 9 74.51 23.6561 30.0461 36.9961 44.5511 52.4911 8.804 12.78 13.9 15.11 15.88 11.32 0.0005
A 10 43.05 27.2373 34.5373 42.7023 51.6873 61.0373 9.0739 14.6 16.33 17.97 18.7 16.621 0.123
B 1 37.8 6.3668 7.6503 9.0298 10.5533 12.2488 2.7893 2.567 2.759 3.047 3.391 3.807 0.3038
B 2 35.4 7.6967 8.9367 10.2797 11.7422 13.3282 3.238 2.48 2.686 2.925 3.172 3.46 0.3769
B 3 29.3 8.6812 10.0817 11.5857 13.2137 14.9497 3.2122 2.801 3.008 3.256 3.472 3.778 0.205
B 4 30.9 8.999 10.3945 11.8675 13.4585 15.18 3.2352 2.791 2.946 3.182 3.443 3.753 0.1478
B 5 29.75 9.5064 10.9409 12.4669 14.1014 15.9644 3.5381 2.869 3.052 3.269 3.726 3.942 0.0644
B 6 39.3 10.1293 11.6418 13.2513 14.9848 16.8748 3.4308 3.025 3.219 3.467 3.78 3.954 -0.0636
B 7 43.3 12.2344 13.8109 15.4949 17.3094 19.3754 3.7561 3.153 3.368 3.629 4.132 4.534 -0.0915
B 8 52 13.2885 14.78 16.4125 18.216 20.2585 4.2523 2.983 3.265 3.607 NA NA -0.0553
B 9 41.6 13.9864 15.6949 17.5539 19.5994 21.7379 4.7464 3.417 3.718 4.091 4.277 4.836 0.0005
B 10 50.2 14.043 16.067 18.2475 20.607 23.178 4.2598 4.048 4.361 4.719 5.142 5.474 0.123
C 1 24.19 1.4639 3.5839 4.8589 6.6539 9.1189 0.3939 4.24 2.55 3.59 4.93 5.63 0.3038
C 2 19.52 2.3384 3.9534 5.7434 7.9884 10.8534 1.3403 3.23 3.58 4.49 5.73 6.51 0.3769
C 3 7.96 2.626 3.711 5.381 7.581 10.846 -1.0251 2.17 3.34 4.4 6.53 5.9 0.205
C 4 15.72 2.7283 4.0083 5.8033 8.0983 11.6433 0.3156 2.56 3.59 4.59 7.09 7.28 0.1478
C 5 11.81 2.3651 3.9451 5.7001 7.9701 10.8651 0.049 3.16 3.51 4.54 5.79 6.39 0.0644
C 6 16.9 2.8771 4.4821 6.5371 9.0671 11.6921 0.964 3.21 4.11 5.06 5.25 7.21 -0.0636
C 7 18.84 3.532 5.792 8.587 12.027 16.417 1.1485 4.52 5.59 6.88 NA NA -0.0915
C 8 21.66 4.5199 7.5199 11.0849 15.4049 20.4099 2.0169 6 7.13 8.64 10.01 11.18 -0.0553
C 9 12.75 4.379 7.449 11.774 16.944 23.244 2.0582 6.14 8.65 10.34 12.6 13.89 0.0005
C 10 24 4.6028 8.4778 13.4228 19.5128 26.4078 0.9412 7.75 9.89 12.18 13.79 16.74 0.123
;
/* Define function XRATE */
proc fcmp outlib=work.funcs.test;
function rs(w, y, B, B1, B2, B3, B4, F1, F2, F3, F4, F5, g, x);
select(w);
when(1) z=-1-exp(-x*1e4); /* for a solution < -1 */
when(2) z=-1+exp(-x*1e4); /* for a solution > -1 */
when(3) z=g+sign((F5-g*B4)*(1+g))*exp(-x*1e4); /* for a solution near g */
end;
p = (F1/(g+z)+F2/(g+z)**2+F3/(g+z)**3+F4/(g+z)**4+F5/(g+z)**5)-(z+g)*(B/(1+g+z)+B1/(1+g+z)**2+B2/(1+g+z)**3+B3/(1+g+z)**4+B4/(1+g+z)**5)+(F5-(g+z)*B4)*(1+g)/(z*(1+g+z)**5);
return(p);
endsub;
function xrate(w, y, p, B, B1, B2, B3, B4, F1, F2, F3, F4, F5, g);
array solvopts[5] initial abconv relconv maxiter status (0.001 0.001 1.0e-6 100);
t=exp(-solve("RS", solvopts, p, w, y, B, B1, B2, B3, B4, F1, F2, F3, F4, F5, g, .)*1e4);
select(w);
when(1) r=-1-t; /* solution < -1 */
when(2) r=-1+t; /* solution > -1 */
when(3) r=g+sign((F5-g*B4)*(1+g))*t; /* solution near g */
end;
return(r);
endfunc;
run;
/* Make function available */
options cmplib=work.funcs;
/* Apply function XRATE to dataset TEST */
data want;
array _B[99] _temporary_;
array _B1[99] _temporary_;
array _B2[99] _temporary_;
array _B3[99] _temporary_;
array _B4[99] _temporary_;
array _P[99] _temporary_;
array _F1[99] _temporary_;
array _F2[99] _temporary_;
array _F3[99] _temporary_;
array _F4[99] _temporary_;
array _F5[99] _temporary_;
do until(last.stock);
set test;
by stock;
r1=xrate(1,y,_P,_B,_B1,_B2,_B3,_B4,_F1,_F2,_F3,_F4,_F5,g);
r2=xrate(2,y,_P,_B,_B1,_B2,_B3,_B4,_F1,_F2,_F3,_F4,_F5,g);
r3=xrate(3,y,_P,_B,_B1,_B2,_B3,_B4,_F1,_F2,_F3,_F4,_F5,g);
output;
end;
Thanks for providing your code.
The new, wide data structure has the array elements B1, B2, ... and F1, F2, ... in one observation, which simplifies the code: We don't need the do until(last.stock) loop any longer. Each observation populates the parameters of the nonlinear equation. There is no need to define individual arrays _B1, _B2, etc.
I had a new idea: It might help to multiply the equation by (r-g)(1+r)**5. Thus we'd get rid of the poles and we'd have to deal only with a polynomial equation in r or perhaps in u:=1+r (of degree 6 in general, with 12 parameters*). This should be numerically more robust. Of course, r=g and r=-1 are still inadmissible.
* Your equation
has parameters B, B1, ..., B4, F1, ..., F5, p and g. Your sample dataset TEST contains a 13th parameter, F (with values that we saw for F1, F2, ... in the equation of the previous thread). Where does this come into play, or can we ignore it?
Thanks for the clarification.
I used the idea mentioned in my previous post and turned the task into finding real zeros of a polynomial in r (of grade 6):
k(r) = c6*r**6 + c5*r**5 + c4*r**4 + c3*r**3 + c2*r**2 + c1*r + c0 = 0
Most of the work went into determining the coefficients c0, ..., c6 as functions of the original 12 parameters. (This is not necessary for computing values of k(r), but it helps to understand the behavior of k(r) -- see below -- and can be useful for various techniques of numerical mathematics that could be applied to it.) To this end I subtracted the right-hand side of your equation from p, multiplied the resulting equation by (r-g)*(1+r)**5 and then collected the terms for r**6, r**5 and so on.
The graphs of several cases (taken from your sample dataset TEST) showed a common pattern: one solution in the interval (-1, 0) and one in the interval (0, 1).
More generally, the constant term of the polynomial,
c0 = g*(B+F1+F2+F3+F4-p)-F5
is negative for each of the 28 observations in TEST with complete sets of parameters (two have "NA" for F4 and F5) and the leading coefficient
c6 = p
is positive throughout your sample data.
That is, k(0)<0 while k(r) goes to infinity for r → ∞ and r → −∞. This implies the existence of at least one positive and at least one negative solution. These are also solutions of the original problem, unless they happen to equal −1 or g (the poles), which should be very unlikely. The definition of function XRATE2 below checks this anyway.
I used −1 (which is a valid argument for the polynomial k, although it's a pole of the original function) and 1 as initial values to find the anticipated solutions around zero.
Here is the code:
/* Define function XRATE2 */
proc fcmp outlib=work.funcs.test;
function k(r, p, B[*], F[*], g);
c6= p;
c5= p*(5-g)-F[1]/2;
c4= p*(10-5*g)-B[1]+B[2]+(g-5)*F[1]/2-F[2]/2;
c3= 10*p*(1-g)+(g-3)*B[1]+(2-g)*B[2]+B[3]+((5*g-9)*F[1]+(g-4)*F[2]-F[3])/2;
c2= p*(5-10*g)+3*(g-1)*B[1]+(1-2*g)*B[2]+(1-g)*B[3]+B[4]+((9*g-7)*F[1]+(4*g-5)*F[2]+(g-3)*F[3]-F[4])/2;
c1= p*(1-5*g)+(3*g-1)*B[1]-g*(B[2]+B[3]+B[4])-(1+g)*B[5]+((7*g-2)*F[1]+(5*g-2)*F[2]+(3*g-2)*F[3]+(g-2)*F[4])/2-F[5];
c0= -p*g+g*B[1]+g*(F[1]+F[2]+F[3]+F[4])-F[5];
z = c6*r**6 + c5*r**5 + c4*r**4 + c3*r**3 + c2*r**2 + c1*r + c0;
return(z);
endfunc;
function xrate2(i, p, B[*], F[*], g); /* first argument is the initial value */
array solvopts[1] i;
r=solve("k", solvopts, 0, ., p, B, F, g);
if min(abs(1+r),abs(r-g))<1e-4 then r=.; /* avoid poles of the original function */
return(r);
endfunc;
run;
/* Make function available */
options cmplib=work.funcs;
/* Apply function XRATE2 to dataset TEST */
data want(rename=(B0=B F0=F));
set test(rename=(B=B0 F=F0));
array B[*] B0-B4;
array F[*] F1-F5;
r1=xrate2(-1, p, B, F, g); /* solution found with initial value -1 */
r2=xrate2( 1, p, B, F, g); /* solution found with initial value 1 */
run;
In each of the 28 complete cases one solution r1 in the interval (−1, 0) and one solution r2 in the interval (0, 1) was found.
I checked all of them using the original rational function definition:
proc fcmp outlib=work.funcs.test;
function h(r, p, B[*], F[*], g);
u=1+r;
s1=p-B[1];
do t=1 to 5;
s1 = s1 - F[t]/u**t;
end;
s2=B[1]/u;
do t=1 to 4;
s2 = s2 + (B[t]+F[t]/2)/u**(t+1);
end;
z = s1 + r*s2 - (F[5]+r*B[5])*(1+g)/((r-g)*u**5);
return(z);
endfunc;
run;
data check;
set want(rename=(B=B0 F=F0));
array B[*] B0-B4;
array F[*] F1-F5;
diff1=h(r1, p, B, F, g);
diff2=h(r2, p, B, F, g);
run;
proc means data=check;
var diff:;
run;
The absolute values of the differences diff1 and diff2 (expected to be close to zero) were <1.5E-8 for your sample data. I'd recommend such a check for all solutions you will find with your real data.
Thank you very much Reinhard,
I will apply the code on other equations, so trying to understand every bit of the procedure here. I got two questions.
First, I understand how c0 (the intercept) is derived. But, having trouble to get C1 to C6. Could you kindly explain it?
Second, in the check procedure, s2 is defined as:
s2 = s2 + (B[t]+F[t]/2)/u**(t+1);
However, the equation has terms of (B[t]+F[t+1]/2)/u**(t+2) rather than (B[t]+F[t]/2)/u**(t+1). Such as,
So, I made an attempt to modify by introducing s3, and separate the above terms into two parts (B[t]/u**(t+2), and (F[t]/2)/u**(t+1)) the code becomes:
proc fcmp outlib=work.funcs.test;
function h(r, p, B[*], F[*], g);
u=1+r;
s1=p-B[1];
do t=1 to 5;
s1 = s1 - F[t]/u**t;
end;
s2=B[1]/u;
do t=1 to 3;
s2 = s2 + B[t]/u**(t+2);
end;
s3 = B[1]/u**2;
do t=1 to 4;
s3 = s3 + (F[t]/2)/u**(t+1);
end;
z = s1 + r*s2 + r*s3 - (F[5]+r*B[4])*(1+g)/((r-g)*u**5);
return(z);
endfunc;
run;
Is my understanding correct? The diff1 and diff2 values become much bigger after this modification.
Thank you.
You're welcome.
@MAC1430 wrote:
First, I understand how c0 (the intercept) is derived. But, having trouble to get C1 to C6. Could you kindly explain it?
After the multiplication by (r-g)*(1+r)**5 all denominators containing r cancel out, but we obtain several terms containing powers of (1+r) and the factor (r-g), now in the "numerator." Multiplying all these out (e.g., (1+r)**4 = 1+4*r+6*r**2+4*r**3+r**4) and aggregating the coefficients of r**j (j=0, ..., 6) we end up with a polynomial in r of grade 6. I did this calculation by hand, but computer algebra software would have helped to obtain the expressions for c0, ..., c6 more easily.
For example, the term F2/(1+r)**2 in the original equation turns into -F2*(r-g)*(1+r)**3 after bringing it to the left side of the equals sign and the multiplication by (r-g)*(1+r)**5. What does this term contribute to, say, c1, i.e., the coefficient of r**1=r? Answer (considering (1+r)**3=1+3*r+3*r**2+r**3): -F2*(-g)*3 - F2*1 = (3*g - 1)*F2. Similarly, the term -r * (B1+0.5*F2)/(1+r)**3 turns into r*(r-g)*(B1+0.5*F2)*(1+r)**2, which contributes -g*(B1+0.5*F2)*1 to c1. Together (focusing on the part containing F2), (3*g - 1)*F2 - g*0.5*F2 = (5*g-2)*F2/2. This is how I obtained the highlighted part of the term for c1:
c1= p*(1-5*g)+(3*g-1)*B[1]-g*(B[2]+B[3]+B[4])-(1+g)*B[5]+((7*g-2)*F[1]+(5*g-2)*F[2]+(3*g-2)*F[3]+(g-2)*F[4])/2-F[5];
@MAC1430 wrote:
Second, in the check procedure, s2 is defined as:
s2 = s2 + (B[t]+F[t]/2)/u**(t+1);
However, the equation has terms of (B[t]+F[t+1]/2)/u**(t+2) rather than (B[t]+F[t]/2)/u**(t+1). Such as,
Note that array B (in the DATA step) is defined as
array B[*] B0-B4;
after renaming variable B from your dataset TEST to B0 in order to avoid a name conflict with the array name. Therefore, B[1]=B0, B[2]=B1, ..., B[5]=B4, whereas the indices of array F are the same as in your formula (F[1]=F1 and so on). I know that this is a bit confusing and maybe I should have named array B differently. Unlike arrays in a DATA step, arrays in PROC FCMP are always 1-based (see the remark about "lower-bound specifications" in the documentation), i.e., a definition like array B[0:4] (to harmonize array indices and original indices) is not allowed.
@MAC1430 wrote:
So, I made an attempt to modify by introducing s3 ...
(...)
Is my understanding correct? The diff1 and diff2 values become much bigger after this modification.
Thank you.
No, it isn't, as explained above, and the substantially increased diff1 and diff2 values demonstrate this. But that's what the check is for: detecting errors.
Thank you very much for the explanation Reinhard.
I applied the procedure on my real dataset, got a mix of positive and negative results in both r1 and r2, some missing values as well. I guess it has something to do with the selection of initial value? Any suggestions how to select an initial value that returns the more robust results? (The unknown r I’m looking for is a rate of return, which should be positive numbers and not exceed 1 or 100%.)
Thank you again for your help, it's really appreciated.
Good to know that only solutions in the interval (0, 1) are of interest. Normally, the algorithm used by the SOLVE function should converge to a particular solution if the initial value is close enough to that solution. Hence, you could use initial values like, for example, 0.01, 0.02, ..., 0.99 in a loop until a solution is found. Or enforce a positive solution by using a substitution similar to what I did in the solution of the earlier thread. Another option would be to switch from SOLVE to the bisection method mentioned in that other post. Can you share an example where you get a negative solution in r2, i.e., with the initial value 1? Then I can have a look which of these options is most suitable.
Thank you very much Reinhard, I selected some sotcks in the attached.
Thanks a lot for the sample data. I checked several cases "manually" and found different scenarios:
I've implemented a search for initial values checking the sign of the polynomial function k at r = 0, 0.001, 0.002, ..., 0.999, 1. The midpoint of an interval with a sign change between its boundaries is selected as an initial value (taking into account even the unlikely case that an interval boundary happens to be a root of k). It turned out that, at least for the 483−67=416 usable sample data observations, the initial values found by this algorithm are close enough to the roots of k in (0, 1), if any, so that the SOLVE function finds the roots accurately. The same roots were found with different subdivisions of [0, 1], e.g., into 10000 subintervals (instead of 1000) or only 3 (!). The interval boundaries enclosing a sign change could also be used as starting values for the bisection method. This would be an alternative approach which could replace the SOLVE function if it didn't work in "extreme" cases.
But I've checked all 171 solutions found for the 416 complete observations (163 with a single solution in (0, 1), 4 with two distinct solutions in (0, 1), 249 with no solutions in (0, 1)) with the code used in the earlier post (i.e., using the rational function h rather than the polynomial), except that no renaming of variables was necessary. These checks resulted in a maximum absolute value of h(r) of 6.08E-6. I took a closer look at the case where this relatively "large" maximum occurred (stock='CH0008742519' & year=1998) and found that
The solutions found with the above approach also confirmed those that I had investigated manually.
Here is the code:
/* Define function XRATE2 */
proc fcmp outlib=work.funcs.test;
function k(r, p, B[*], F[*], g);
c6= p;
c5= p*(5-g)-F[1]/2;
c4= p*(10-5*g)-B[1]+B[2]+(g-5)*F[1]/2-F[2]/2;
c3= 10*p*(1-g)+(g-3)*B[1]+(2-g)*B[2]+B[3]+((5*g-9)*F[1]+(g-4)*F[2]-F[3])/2;
c2= p*(5-10*g)+3*(g-1)*B[1]+(1-2*g)*B[2]+(1-g)*B[3]+B[4]+((9*g-7)*F[1]+(4*g-5)*F[2]+(g-3)*F[3]-F[4])/2;
c1= p*(1-5*g)+(3*g-1)*B[1]-g*(B[2]+B[3]+B[4])-(1+g)*B[5]+((7*g-2)*F[1]+(5*g-2)*F[2]+(3*g-2)*F[3]+(g-2)*F[4])/2-F[5];
c0= -p*g+g*B[1]+g*(F[1]+F[2]+F[3]+F[4])-F[5];
z = c6*r**6 + c5*r**5 + c4*r**4 + c3*r**3 + c2*r**2 + c1*r + c0;
return(z);
endfunc;
function xrate2(i, p, B[*], F[*], g); /* first argument is the initial value */
array solvopts[1] i;
r=solve("k", solvopts, 0, ., p, B, F, g);
if min(abs(1+r),abs(r-g))<1e-4 then r=.; /* avoid poles of the original function */
return(r);
endfunc;
run;
/* Make function available */
options cmplib=work.funcs;
/* Find initial values for solutions in (0, 1) */
%let nsubintv=1000; /* number of subintervals of [0, 1] to search */
data init(drop=_:);
set have;
if nmiss(P, g, of F1-F5, of B0-B4) then output; /* Solutions will be missing for incomplete data. */
else do;
array init[6]; /* Up to six roots are possible theoretically. */
array B[*] B0-B4;
array F[*] F1-F5;
do _i=0 to &nsubintv;
_x=_i/&nsubintv;
_z=k(_x, p, B, F, g);
_s=sign(_z);
_prev_s=lag(_s);
if _s~=0 then do;
if _i & _s~=_prev_s~=0 then do;
_j=sum(_j,1);
init[_j]=(_i-0.5)/&nsubintv; /* select midpoint of subinterval with sign change between boundaries */
end;
end;
else do; /* This case (solution on subinterval boundary) is very unlikely. */
_j=sum(_j,1);
init[_j]=_x;
end;
end;
output;
end;
run;
/* Apply function XRATE2 to dataset INIT */
data want(drop=_j);
set init;
array B[*] B0-B4;
array F[*] F1-F5;
array init[6];
array r[6];
do _j=1 to dim(r) while(init[_j]>.);
r[_j]=xrate2(init[_j], p, B, F, g); /* solution found with initial value init[_j] */
end;
run;
The PROC FCMP step and the OPTIONS statement are exactly the same as in my earlier post. I repeated them because your attachment sample.sas contained something different. For your sample data the dimension of the arrays init and r in the DATA steps could be reduced from 6 to 2 (remaining variables are always missing), but theoretically a polynomial of grade 6 can have up to 6 real roots, so it's better to allow for that many solutions and initial values. You may want to drop variables init1-init6. I've left them in dataset WANT just to check the differences (always within ±0.0005 for &nsubintv=1000) between them and the corresponding solutions.
Are you ready for the spotlight? We're accepting content ideas for SAS Innovate 2025 to be held May 6-9 in Orlando, FL. The call is open until September 25. Read more here about why you should contribute and what is in it for you!
Learn how use the CAT functions in SAS to join values from multiple variables into a single value.
Find more tutorials on the SAS Users YouTube channel.