Solved: Re: FCMP - solving equation

MAC1430 · Posted 09-15-2021 02:50 AM

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).

FreelanceReinh · Posted 09-21-2021 03:00 PM

Thanks a lot for the sample data. I checked several cases "manually" and found different scenarios:

no real (as opposed to complex) solution at all (e.g., stock='BE0003826436', year=2012, where the polynomial function is everywhere positive while the original function skips the x-axis at two poles!)
only two negative solutions (not always found by SOLVE with initial values −1 or 1)
one positive solution >1 and one negative solution that SOLVE found with both initial values
one negative solution and one in (0, 1) that SOLVE found with both initial values
one negative solution and one in (0, 1), both found by SOLVE
two positive solutions: one in (0, 1) and one >1, but sometimes only the former was found by SOLVE (with both initial values, −1 and 1)
two solutions in (0, 1), both found by SOLVE
missing solutions, of course, where one or more parameters used in the formula are missing ("NA") -- this applies to 67 of the 483 observations in the sample data.

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 root determined by SOLVE differs from the "best possible" approximation (of the exact solution) by only <5E-10, which is well within expectations.
The value of h(r) could be improved to <8E-13 by changing r by the above tiny amount.
The first derivative of h at the solution r=0.8968... (which is the only solution in (0, 1) here) is >14,000 (!). This explains the observed "sensitivity" to minor changes of r.

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.

View solution in original post

PeterClemmensen · Posted 09-15-2021 02:55 AM

Did you try to modify @FreelanceReinhs solution?

The DATA to DATA Step Macro
Blog: SASnrd

MAC1430 · Posted 09-15-2021 03:03 AM

yes, I tried and still trying, got error "ERROR 71-185: The xrate function call does not have enough arguments".

FreelanceReinh · Posted 09-15-2021 05:46 AM

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.

MAC1430 · Posted 09-16-2021 12:05 AM

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;

FreelanceReinh · Posted 09-16-2021 04:05 AM

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?

MAC1430 · Posted 09-16-2021 09:15 PM

Thank you very much Reinhard, F was used to calculate F1 to F5, yes, we can ignore F for this equation.

FreelanceReinh · Posted 09-17-2021 10:13 AM

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.

MAC1430 · Posted 09-19-2021 06:42 AM

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.

FreelanceReinh · Posted 09-19-2021 09:53 AM

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.

MAC1430 · Posted 09-20-2021 08:00 AM

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.

FreelanceReinh · Posted 09-20-2021 09:28 AM

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.

MAC1430 · Posted 09-21-2021 01:02 AM

Thank you very much Reinhard, I selected some sotcks in the attached.

FreelanceReinh · Posted 09-21-2021 03:00 PM

Thanks a lot for the sample data. I checked several cases "manually" and found different scenarios:

no real (as opposed to complex) solution at all (e.g., stock='BE0003826436', year=2012, where the polynomial function is everywhere positive while the original function skips the x-axis at two poles!)
only two negative solutions (not always found by SOLVE with initial values −1 or 1)
one positive solution >1 and one negative solution that SOLVE found with both initial values
one negative solution and one in (0, 1) that SOLVE found with both initial values
one negative solution and one in (0, 1), both found by SOLVE
two positive solutions: one in (0, 1) and one >1, but sometimes only the former was found by SOLVE (with both initial values, −1 and 1)
two solutions in (0, 1), both found by SOLVE
missing solutions, of course, where one or more parameters used in the formula are missing ("NA") -- this applies to 67 of the 483 observations in the sample data.

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 root determined by SOLVE differs from the "best possible" approximation (of the exact solution) by only <5E-10, which is well within expectations.
The value of h(r) could be improved to <8E-13 by changing r by the above tiny amount.
The first derivative of h at the solution r=0.8968... (which is the only solution in (0, 1) here) is >14,000 (!). This explains the observed "sensitivity" to minor changes of r.

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.

MAC1430 · Posted 09-22-2021 04:17 AM

Thank you very much for everything, you are amazing.

Ready to join fellow brilliant minds for the SAS Hackathon?

Classroom Training Available!