Best 5 segment piecewise linear fit

Question

I'm trying to get a piecewise linear best-fit for the closing price of one of the stocks I'm interested in. The logic seems ok, and the workflow works for a straight-line (ie 2 pts, ie 1 segment) regression (actually NMinimize)... but if I increase the the number of variables to solve for, it breaks and complains about "abscissa"... the code is below with the comments. Also here is the data file . Can you help me figure out whats wrong? (also here is the nb file in case you need it).

In[120]:= (d = 
   Import["data_out.txt", "CSV"]);

In[121]:= d = 
 d[[4533 ;; 4922, 
   5]]; (*get close price data for 14-may from the dataset*)

In[122]:= d = 
 Flatten[{Table[i, {i, 390}], 
   d}, {2}]; (*add an index for the price ie 1-390 data pts*)

In[123]:= (* create objective function to minimize:
1. piecewise linear interpoltion function takes a set of points "p" 
2. applies function to the index
3. subtracts the close price 
4. squares the diffrenence
5. sum
6. root *)

In[124]:= 
e[p_] := Total[(Interpolation[p, InterpolationOrder -> 1]@
       d[[All, 1]] - d[[All, 2]])^2]^0.5

In[125]:= (*this is the solution i got using excel solver*)

In[126]:= 
excelsolution = {{1, 32.69967765}, {28.16280834, 
    31.37817608}, {108.0001043, 32.75429029}, {135.5658831, 
    31.7584233}, {299.8762066, 32.76192525}, {390, 32.88427106}};

In[127]:= e[excelsolution]

Out[127]= 2.13146

In[128]:= ListPlot[{d, excelsolution}, Joined -> {False, True}, 
 PlotMarkers -> {{Automatic, Tiny}, {Automatic, Small}}] (* run to see it *)

Out[128]= (*graphics pasted above*)

In[129]:= (*this works*)

In[130]:= NMinimize[
 {
  e[{{1, y0}, {390, y390}}],
  31.5 <= y0 <= 33 && 31.5 <= y390 <= 33
  },
 {y0, y390},
 Method -> "DifferentialEvolution"
 ]

Out[130]= {5.51241, {y0 -> 31.7968, y390 -> 32.8737}}

In[131]:= (*but this doesn't work when i increase the number of \
arguments*)

In[132]:= NMinimize[
 {
  e[{{1, y0}, {x1, y1}, {390, y390}}], 
  31.5` <= y0 <= 33 && 1 <= x1 <= 390 && 31.5` <= y1 <= 33 && 
   31.5` <= y390 <= 33
  },
 {y0, x1, y1, y390},
 Method -> "DifferentialEvolution"
 ]

During evaluation of In[132]:= Interpolation::indat: Data point {x1,y1} contains abscissa x1, which is not a real number.

During evaluation of In[132]:= Interpolation::indat: Data point {x1,y1} contains abscissa x1, which is not a real number.

During evaluation of In[132]:= Interpolation::indat: Data point {x1,y1} contains abscissa x1, which is not a real number.

During evaluation of In[132]:= General::stop: Further output of Interpolation::indat will be suppressed during this calculation.

During evaluation of In[132]:= NMinimize::nnum: The function value {10.8271,10.8107,10.7944,10.7781,10.7618,10.7456,10.7294,10.7133,<<35>>,10.157,10.1423,10.1276,10.113,10.0984,10.0839,10.0694,<<340>>} is not a number at {x1,y0,y1,y390} = {285.476,32.6896,32.3246,32.9648}.

Answer 1

Here is an approach based on generating a Line object, determining the total distance from all points in your dataset to the line, and minimizing that distance. In the following, price is obtained from your data, similar to what you described in your code:

price = Import["data_out.txt", "CSV"][[4533 ;; 4922, 5]];
price = Transpose@{Range[Length[price]], price}

First, let's define a target function to minimize:

ClearAll[model]

model[data_, positionlist_] /;
  (Max[positionlist] >= Length[data] || Min[positionlist] < 1) := 10.^10

model[data_, positionlist_?(VectorQ[#, NumericQ] &)] := Module[{rdf},
   rdf = RegionDistance[Line@data[[Join[{1}, Sort@positionlist, {-1}] ]] ];
   Total@rdf[data]
 ]

The minimization is carried out using NArgMin using the Simulated Annealing algorithm with a very high perturbation scale setting to explore a wider swath of the parameter space. Such a high value of the perturbation scale generates quite a few "unacceptable" values (e.g. negative, or larger than the size of the dataset); rather than introducing a constraint, which slows down the minimization significantly, I introduced "guard definitions" in the target model function, as shown above, that return an extremely high value.

min = Sort@Round@
   NArgMin[
     model[price, Round@{a, b, c, d}],
     {a, b, c, d},
     Method -> {
       "SimulatedAnnealing",
       "PerturbationScale" -> 100
     }
   ]

Here is the resulting line overlaid on the points:

ListPlot[
  price,
  Epilog -> {
    PointSize[0.02],
    Through[{Point, Line}@price[[{1, Sequence @@ min, -1}]]]
  }
]

---

Although this approach does not require a manual choice of starting conditions as required, I want to point that many similar minima are present, so even very small changes in any of the parameters may provide rather significantly different results. Unfortunately, this seems to be a feature of the problem, rather than of the minimization process.

Answer 2

The reason you are receiving the "Data point {x1,y1} contains abscissa x1, which is not a real number" error is that Interpolation expects a numerical value for the abscissa.

In the first example for the end points you used 1 and 390. This worked fine.

In the second example you used, in addition, x1.

That was what Interpolation was complaining about.

It may be overdoing it but functions were defined so that all arguments were required to be numeric.

Initialization

dataOut = Import["data_out.txt", "CSV"];
d = dataOut[[4533 ;; 4922, 5]];

A function was created using the Excel results.

intFunExcel = Interpolation[excelsolution, InterpolationOrder -> 1]

and then the data and Excel function were plotted.

Show[
 ListPlot[
  d,
  PlotStyle -> Red,
  PlotMarkers -> {Automatic, Tiny}
  ],
 Plot[
  intFunExcel[x],
  {x, 1, 390},
  PlotStyle -> Black
  ]
 ]

NMinimize - numeric functions

This was broken into two parts, first the function that will return the y value and then secondly the error function.

Given an x value, curveFitIntF will return a y value. It needs the x input value and six data points (the first and last x values are 1 and 390).

curveFitIntF[
  x_?NumericQ,
   y1_?NumericQ,
  {x2_?NumericQ, y2_?NumericQ},
  {x3_?NumericQ, y3_?NumericQ},
  {x4_?NumericQ, y4_?NumericQ},
  {x5_?NumericQ, y5_?NumericQ},
  y390_?NumericQ
  ] := Interpolation[
   {
    { 1.0, y1}, 
    {x2, y2},
    {x3, y3},
    {x4, y4},
    {x5, y5},
    {390.0, y390}
    },
   InterpolationOrder -> 1][x]

Next define the error function. Note: The error itself rather than square root was used.

error[ y1_?NumericQ,
  {x2_?NumericQ, y2_?NumericQ},
  {x3_?NumericQ, y3_?NumericQ},
  {x4_?NumericQ, y4_?NumericQ},
  {x5_?NumericQ, y5_?NumericQ},
  y390_?NumericQ] := Total[
  Map[
   (curveFitIntF[#[[1]], y1, {x2, y2}, {x3, y3}, {x4, y4}, {x5, y5}, 
        y390] - #[[2]])^2 &,
   d]
  ]

Next, run it through NMinimize.

Note: NMinimize really needs some reasonable starting points for the parameters. This is quite important.

I eyeballed them from the graph and used +/- 10 for the x value and (much smaller) +/- 0.1 for the y value.

seg5 = NMinimize[
  {
   error[
     y1,
    {x2, y2},
    {x3, y3},
    {x4, y4},
    {x5, y5},
    y390
    ],
   1.0 < x2 < x3 < x4 < x5 < 390.0 &&
    31 < y1 < 33 &&
    31 < y2 < 33 &&
    31 < y3 < 33 &&
    31 < y4 < 33 &&
    31 < y5 < 33 &&
    31 < y390 < 33
   },
  {
   {y1, 32.6, 32.8},
   {x2, 25.0, 35.0},
   {y2, 31.3, 31.5},
   {x3, 90.0, 110.0},
   {y3, 32.7, 32.9},
   {x4, 110, 130},
   {y4, 31.6, 31.8},
   {x5, 290, 310},
   {y5, 32.6, 32.8},
   {y390, 32.7, 32.9}
   },
  Method -> "DifferentialEvolution"
  ]

The results were:

{4.51966, {y1 -> 32.6719, x2 -> 27.7708, y2 -> 31.3766, x3 -> 107.248,
   y3 -> 32.7467, x4 -> 136.495, y4 -> 31.7585, x5 -> 299.32, 
  y5 -> 32.7626, y390 -> 32.8829}}

This produced an error very slightly smaller than the excel results.

A function is defined to use the seg5 results.

intFun = Interpolation[
   {
    { 1.0, y1}, 
    {x2, y2},
    {x3, y3},
    {x4, y4},
    {x5, y5},
    {390.0, y390}
    },
   InterpolationOrder -> 1] /. seg5[[2]]

One can ignore the warning message.

Now plot it:

Show[
 ListPlot[
  d,
  PlotStyle -> Red,
  PlotMarkers -> {Automatic, Tiny}
  ],
 Plot[
  intFun[x],
  {x, 1, 390},
  PlotStyle -> Black
  ]
 ]