Solving an equation , non linear

[Thomas_mp]

Good morning,

I am trying to solve an equation of the type   a1 = c1+  (b1-X)*c1/(1+X)  + (b2-X)*c2/(1+X)^2 + (b3-X)*c3/(1+X)^3 +… (bt-X)*ct/X*(1/(1+X)^t-1)

Note that the first (t-1) terms are of the same form (bj-Xj)*cj/(1+X)^j-1  , but the last term is different.

I need to find the value of X ( 0<X<1).  In the case of multiple solutions, chose the smaller non-negative solution. 

In my case I have 25 variables. 

• ai
• 12 bi
• 12 ci

The data looks like:

Input   a    b1    b2    b3    b4    b5    b6    b7    b8    b9    b10  b11    b12         c1    c2    c3    c4    c5    c6    c7    c8    c9    c10  c11    c12              

cards

84.3  0.33  0.31  0.29  0.28  0.26  0.25  0.24  0.22  0.21  0.19  0.18  0.17        11.1  13.4  16.1  19.1  22.4  26.2  30.4  34.9  39.8  45.1  50.7  56.5

91.5  0.332  0.310  0.291  0.277  0.263  0.249  0.235  0.221  0.207  0.194  0.180  0.166  11.076  13.422  16.074  19.061  22.433  26.201  30.368  34.927  39.859  45.132  50.700  56.502

I can, by trial and error, find the value of X, this is: solve the equation  in Excel

84.3=  11.1 +(0.33-X)*11.1/(1+X) +  (0.31-X)*13.4/(1+X)^2 + (0.29-X)*16.1/(1+X)^3+ (0.28-X)*19.1/(1+X)^4 + .............  + (0.18-X)*50.7/(1+X)^11 +  (0.17-X)*56.5/X*[1/(1+X)^11]

In this case the solution is X=0.068, this X produces a value for  a=84.3.

Doing the same, the solution for the second row of data (which starts with a value of a=91.5) is X = 0.0638

But I have several thousands of observations …  so I do not how to do this in Excel by trial and error (function IRR does not work in excel because X is both in the numerator and in the denominator).

I am attaching the excel file with the data and the solution. 

Thank you !!!

Tomas

[RobPratt · Accepted solution]

Here's one way to do it with PROC OPTMODEL by using a DO loop to call the solver once per observation:

%let n = 12;
data indata;
   input a b1-b&n c1-c&n @@;
   datalines;
84.3 0.33  0.31  0.29  0.28  0.26  0.25  0.24  0.22  0.21  0.19  0.18  0.17  
     11.1  13.4  16.1  19.1  22.4  26.2  30.4  34.9  39.8  45.1  50.7  56.5
91.5 0.332 0.310 0.291 0.277 0.263 0.249 0.235 0.221 0.207 0.194 0.180 0.166 
     11.076 13.422 16.074 19.061 22.433 26.201 30.368 34.927 39.859 45.132 50.700 56.502
;

proc optmodel printlevel=0;
   num n = &n;
   set OBS;
   num a {OBS};
   num b {OBS, 1..n};
   num c {OBS, 1..n};
   read data indata into OBS=[_N_]
      a {j in 1..n} <b[_N_,j]=col('b'||j)> {j in 1..n} <c[_N_,j]=col('c'||j)>;

   var X >= 0 <= 1;
   min Z = X;
   num obsThis, aThis, bThis {1..n}, cThis {1..n};
   con IRR:
      aThis = cThis[1] + sum {j in 1..n-1} (bThis[j]-X)*cThis[j]/(1+X)^j + (bThis[n]-X)*cThis[n]/X/(1+X)^(n-1);

   num Xsol {OBS};
   do obsThis = OBS;
      put obsThis=;
      aThis = a[obsThis];
      for {j in 1..n} do;
         bThis[j] = b[obsThis,j];
         cThis[j] = c[obsThis,j];
      end;
      solve;
      Xsol[obsThis] = X.sol;
   end;
   create data outdata from [i] X=Xsol;
quit;

You can also use a COFOR loop to solve these independent problems concurrently.  Just replace the DO statement with these two statements:

   cofor {i in OBS} do;
      obsThis = i;

In SAS Viya 3.5, it is even simpler (no explicit looping) by using the groupBy option in the runOptmodel action:

data sascas1.indata;
   set indata;
   i = _N_;
run;

proc cas noqueue;
   source pgm;
      put _BY_LINE_;

      num n = &n;
      num a;
      num b {1..n};
      num c {1..n};
      read data indata into
         a {j in 1..n} <b[j]=col('b'||j)> {j in 1..n} <c[j]=col('c'||j)>;

      var X >= 0 <= 1;
      min Z = X;
      con IRR:
         a = c[1] + sum {j in 1..n-1} (b[j]-X)*c[j]/(1+X)^j + (b[n]-X)*c[n]/X/(1+X)^(n-1);

      solve;
      create data outdata from X;
   endsource;
   
   action optimization.runOptmodel / code=pgm printlevel=0 groupBy={"i"};
   run;
quit;

[Thomas_mp]

Thank you !!!

[Thomas_mp]

Robert, Your code using proc optmodel works. Thank you. I am trying to merge the solutions for X (in outdata) with my original sample (data - indata). I can create a counting variable n_ +1 in indata.
How can I keep this variable n_ in outdata so I can merge outdata and indata by n_ ?

Thank you again
Tomas

[RobPratt]

No need to merge.  Just replace the READ DATA statement with this:

read data indata into OBS=[n_]
   a {j in 1..n} <b[n_,j]=col('b'||j)> {j in 1..n} <c[n_,j]=col('c'||j)>;

And then replace the CREATE DATA statement with this:

create data outdata from [n_] X=Xsol
   a {j in 1..n} <col('b'||j)=b[n_,j]> {j in 1..n} <col('c'||j)=c[n_,j]>;

[Thomas_mp]

Thank you again !!!