nl using a piecewise function

[Viridiana Corona]

Dear stata community:

I am trying to fit several data sets with three distinct regions. Essentially, I want to fit the coefficients (a1,a2,a3,b1,b2) of the following piecewise function, where I also want to determine the best kinks location (d,e):



The following Stata page explains how to find the coefficients and kinks for a function in two parts: (https://stats.idre.ucla.edu/stata/fa...se-regression/), but I have problems setting up the function with 3 regions. My problem is how do I define the x's in my nl function to consider the kinks? Here is what I have tried:

nl (y = ({a1}*(x^({b1}-1))-{c})*(x < {d})+ ///
        (({a1}*({d}^({b1}-1))-{d}) + {a2}*(x-{d}))*({d} <=x < {e})+///
        ({a2}*{e} + {a3}*(x-{e})^({b2}-1))*({x >= {e})),///
        initial(a1 0 a2 0 a3 0 b1 1.1 b2 2.1 c 50 d 18000 e 25000)

Also, from what I have read, using nl is a bit complicated because its results depend greatly on the initial values. Do you have any alternative to solve this problem? Maybe linearizing the first and second parts is easier for Stata? However, the parameter c in the first part makes it impossible to split the linearized function.

Thank you very much for all your help!

I attached a sample image and a sample dataset:

* Example generated by -dataex-. To install: ssc install dataex
clear
input float(y x)
   -3000 17305.16
-1888.98 17305.16
 -1888.9 17317.16
   -1489 17317.16
 -1488.9 17342.16
    -889 17342.16
  -888.9 17367.16
  -500.1 17367.16
    -500 17655.16
  -400.1 17656.16
    -400 17675.16
    -339 17676.16
  -338.9 17701.16
  -300.1 17701.16
    -300 17711.16
  -200.1 17712.16
    -200 17732.16
 -199.98 17732.16
  -199.9 17770.16
    -166 17771.16
  -165.9 17771.16
  -150.1 17771.16
    -150 17779.16
    -139 17779.16
  -138.9 17829.16
    -100 17829.16
   -99.9 17829.16
     -51 17830.16
   -50.1 17858.16
     -50 17866.16
   -49.9 17866.16
     -36 17866.16
   -35.9 17891.16
   -30.1 17891.16
     -30 17893.16
     -25 17893.16
   -24.9 17893.16
     -20 17893.16
   -19.9 17895.16
   -10.8 17895.16
   -10.1 17895.16
     -10 17908.16
    -9.9 17912.16
      -6 17912.16
    -5.9 17912.16
    -5.1 17912.16
      -5 17912.16
    -3.1 17912.16
      -3 17912.16
    -2.5 17912.16
    -2.4 17940.16
    -1.5 17940.16
    -1.4 17940.16
    -1.1 17940.16
      -1 17940.16
    -.57 17940.16
     -.5 17970.16
     -.1 17970.16
       0 18035.16
      .1 18255.16
      .2 18255.16
      .9 18255.16
       1 18268.16
     1.9 18268.16
       2 18271.16
     2.1 18293.16
       3 18293.16
       4 18294.16
       5 18294.16
     5.1 18294.16
     5.6 18294.16
     5.7 18294.16
       6 18294.16
     6.1 18294.16
       9 18294.16
     9.2 18294.16
     9.3 18294.16
     9.9 18294.16
      10 18294.16
    10.1 18319.16
      11 18319.16
    11.1 18319.16
    11.9 18320.16
      12 18350.16
    12.1 18350.16
      13 18350.16
      14 18350.16
    14.9 18353.16
      15 18354.16
      16 18355.16
    16.9 18355.16
      17 18360.16
    17.6 18360.16
    17.7 18360.16
      18 18360.16
    18.7 18360.16
    18.8 18395.16
      19 18395.16
      19 18395.16
      20 18437.16
end

[Viridiana Corona]

I think I found a way to define the function using nl. However, now I get the error "starting values invalid or some RHS variables have missing values". Can someone, please, help me to know why? Thank you!

nl (y = ({a1}*(x^({b1}-1))-{c})*(x < {e}) + ({a1}*({e}^({b1}-1))-{c} + {a2}*(x-{e}))*({e}<= x < {d})+ ({a2}*({d}-{e}) + {a3}*((x-{d})^({b3}-1)))*(x >= {d})), initial( a1 0.1 a2 1 a3 0.1 b1 1.1 b3 2.1 e 17000 d 24000)