Programmatic Identification Of SupportResistance Trend Lines With Python

Trend line analysis to identify support and resistance levels was traditionally done by economists by hand drawing lines on charts such as the closing price chart for a particular security. The computerized automation of such a task has widely not been properly implemented in a great deal of libraries out there. The start of such a task requires getting an objective definition of trend lines:

“In finance, a trend line is a bounding line for the price movement of a security. It is formed when a diagonal line can be drawn between a minimum of three or more price pivot points. A line can be drawn between any two points, but it does not qualify as a trend line until tested. Hence the need for the third point, the test.”

The three or more point requirement is where almost all trivial attempts fall short. So we will discuss the technical details behind each aspect of implementing this. This could further be used to identify various chart pattern such as triangular shapes which includes triangles, wedges, pennants and parallel shapes such as flags, double/triple tops/bottoms, head and shoulders, rectangles although all of these often require only 2 points, yet would doubtlessly have a strengthened indication if found with 3 where applicable. This interesting topic will be saved for future discussion.

We will use Python 3 (though any version from 2.7 on should be adequate) for this task as it has good libraries designed to deal with data sets and is quite easy to use. First, the stock data must be fetched. For this series of examples, the prevalent index of the S&P500 (Ticker: ^GSPC) will be used to have practical relevance.

The Python library yfinance can very easily fetch the entire historical data. For example purposes it is convenient to take slightly less than 4 years of data or the last 1000 points. It will return the data as a pandas DataFrame for convenient indexing and other operations:

import numpy as np
import pandas as pd
import yfinance as yf #pip install yfinance
tick = yf.Ticker(‘^GSPC’)
hist = tick.history(period=”max”, rounding=True)
#hist = hist[:’ ‘]
hist = hist[-1000:]
h = hist.Close.tolist()

For reproducible results, the data set was taken as of closing on October 7, 2019 and a comment filtering to that effect is shown.

The first issue is to identify pivot points. Our points are going to be the closing prices at a given time. We can refer to these points in the chart as peaks and troughs or as local maxima and local minima.

All peaks and troughs naive method
There is a naive method for doing this as a pivot point would require the point preceding and succeeding to be both lower or both higher than a current point. A naive method has serious shortcomings however, as if the price remained constant for 2 consecutive days, no peak or trough would be detected. Nonetheless the indexes where this scenario did not occurred can be calculated relatively easily:

minimaIdxs = np.flatnonzero(
hist.Close.rolling(window=3, min_periods=1, center=True).aggregate(
lambda x: len(x) == 3 and x[0] > x[1] and x[2] > x[1])).tolist()
maximaIdxs = np.flatnonzero(
hist.Close.rolling(window=3, min_periods=1, center=True).aggregate(
lambda x: len(x) == 3 and x[0]

All peaks and troughs handling consecutive duplicates method
Certainly, collapsing the places where the price remains constant and the prior code with a different index retrieval would overcome this shortcoming and yield all of them, with a single time in the collapsed flat interval identified. Dropping consecutive duplicates before the above calculation is as simple as: hist.Close.loc[hist.Close.shift(-1) != hist.Close]. However indexes need re-computation as they can have changed:

hs = hist.Close.loc[hist.Close.shift(-1) != hist.Close]
x = hs.rolling(window=3, center=True)
.aggregate(lambda x: x[0] > x[1] and x[2] > x[1])
minimaIdxs = [hist.index.get_loc(y) for y in x[x == 1].index]
x = hs.rolling(window=3, center=True)
.aggregate(lambda x: x[0] maximaIdxs = [hist.index.get_loc(y) for y in x[x == 1].index]

Both of these methods can be done using a list and loop as opposed to using the Pandas rolling window technique and would be significantly faster in fact with that simple implementation change. Although Pandas provides short elegant code snippets, its not always the most efficient to utilize them especially the aggregate function which uses callbacks.

Numerical Differentiation method
However a far more elegant approach to this issue is to use the numerical derivative of the closing price to identify the points. The first derivative is the rate of change or effectively the momentum or velocity of the closing price, while the second derivative indicates the rate of change of the first derivative or its acceleration. Normal differentiation does not apply here as a discrete time series requires discrete numerical analysis tools. There are several advantages including that the numerical derivative will smooth out the data by considering all points in a given range from the point at which the rate of change is to be computed. For example the 5 points stencil method which considers some increment before and ahead, as well as a double increment before and ahead along with the point itself. The findiff library makes this computation easy and accurate even using higher order approximation methods:

from findiff import FinDiff #pip install findiff
dx = 1 #1 day interval
d_dx = FinDiff(0, dx, 1)
d2_dx2 = FinDiff(0, dx, 2)
clarr = np.asarray(hist.Close)
mom = d_dx(clarr)
momacc = d2_dx2(clarr)

With the first and second derivative calculated, a degree of smoothing has effectively occurred giving us the prominent peaks and troughs. They are the places where the first derivative is 0 as no momentum indicates a change of direction has occurred. The second derivative being positive or negative indicates a trough or peak, respectively as upward verses downward acceleration indicates the reversal in direction. However, an exact first derivative of 0 is very unlikely. Instead a value on one side of 0 followed by another is realistically what will occur since the point of 0 derivative occurred in between two days. So based on this, the higher or lower closing price will be chosen between two points where a 0 crossover has occurred:

def get_extrema(isMin):
return [x for x in range(len(mom))
if (momacc[x] > 0 if isMin else momacc[x] (mom[x] == 0 or #slope is (x != len(mom) – 1 and #check next day (mom[x] > 0 and mom[x+1] h[x] >= h[x+1] or mom[x] 0 and h[x] x != 0 and #check prior day (mom[x-1] > 0 and mom[x] h[x-1] mom[x-1] 0 and h[x-1] > h[x])))]
minimaIdxs, maximaIdxs = get_extrema(True), get_extrema(False)

Quite a lengthy discussion, and if all minima and maxima points are desired, then combining the indexes from the naive method could be done as it would capture ones that were smoothed out during the momentum computation. But for the purpose of trend lines, the prominent pivot points are typically ideal and the other ones are mostly noisy or irrelevant. Here is a small plot showing 10 days with their velocity and acceleration values along with the pivot points identified.

Closing Price with Pivot Points, Momentum, AccelerationNow for the mathematically curious who do not want to leave this computation merely up to a library like findiff, I will suggest how to calculate a single point in this data (2019–09–27 closed at 2961.79, its prior was 2977.62 and its succeeding 2976.74). It computes its own coefficients to enable higher order methods.

Most people understand derivative as the change (delta — △) of y over the change of x. For a continuous line this is merely the slope of the line and with subsequent points it is trivial to calculate by taking a difference of y values. But the derivative is actually the limit of △x over △y as △x approaches 0. Taking such a limit on discrete data points requires actually broadening the amount of points being looked at.

This means, there is technically a data leak here though its a trivial and non-concerning one, as future values are considered for past data. It is subtle and not particularly important, though if finding trend lines using very recent data within the time interval window of the numerical derivative, perhaps a different technique mentioned is better. How many days are we talking about? Well it depends on the accuracy, by default only 1 day on each side of the window, except for the left most and right most values which would look 2 days ahead or behind.

So I will give code showing hand calculations of one central point:

import findiff
coeff = findiff.coefficients(deriv=1, acc=1)
print(coeff)

This will yield:

{‘center’: {‘coefficients’: array([-0.5, 0. , 0.5]),
‘offsets’: array([-1, 0, 1])},
‘forward’: {‘coefficients’: array([-1.5, 2. , -0.5]),
‘offsets’: array([0, 1, 2])},
‘backward’: {‘coefficients’: array([ 0.5, -2. , 1.5]),
‘offsets’: array([-2, -1, 0])}}

So findiff uses coefficients of -0.5, 0, 0.5 and time differences of -1, 0 and 1 for all the center points. Of course on the far right and far left, it uses the forward and backward values. Calculating all of the derivatives is straight forward then. The details of calculating what is termed the finite difference coefficients and window sizes is outside the scope here but tables and methods for easily computing them exist.

hist = hist[:’2019–10–07’]
day = # September 27, 2019 per example (-7 or 7th last point)
sum([coeff[‘center’][‘coefficients’][x] *
hist.Close[day + coeff[‘center’][‘offsets’][x]]
for x in range(len(coeff[‘center’][‘coefficients’]))])

This yields as the picture shows:

* -0. =-2977.62/2+2976.74/2

For the acceleration:

coeff=findiff.coefficients(deriv=2, acc=1)
print(coeff)

Yields:

{‘center’: {‘coefficients’: array([ 1., -2., 1.]),
‘offsets’: array([-1, 0, 1])},
‘forward’: {‘coefficients’: array([ 2., -5., 4., -1.]),
‘offsets’: array([0, 1, 2, 3])},
‘backward’: {‘coefficients’: array([-1., 4., -5., 2.]),
‘offsets’: array([-3, -2, -1, 0])}}

And again:

sum([coeff[‘center’][‘coefficients’][x] *
hist.Close[day + coeff[‘center’][‘offsets’][x]]
for x in range(len(coeff[‘center’][‘coefficients’]))])

And we also see the expected value:

* 30. =2977.62–2961.79*2+2976.74

The accuracy can be increased though it looks like the default of 1 is used so its not the 5 point stencil I pointed you too but has its own technique to generate the offsets and coefficients which Wikipedia also goes into some detail about. The 5 point stencil coefficients which have already been divided there can be found with: print(findiff.coefficients(deriv=1, acc=3)). Then it looks 2 days ahead and behind instead of just 1 with [1/12, -2/3, 0, 2/3, -1/12] identical and merely in reverse order as presented. Note that acc=3 can be passed as an additional parameter to the FinDiff constructor for more accuracy.

Having finally gotten the pivot points, selecting 2 of them would be very easy as any pair of the minima or maxima points would constitute a line for a total of n(n-1) total lines which are all perfect fitting lines as basic geometry indicates that 2 unique points is one description of a line. However 3 points could be selected in n(n-1)(n-2) ways and such enumeration is of algorithmic complexity O(n³). These points will not always be on a straight line so there will also be some degree of error. The degree of error will be important since we need to filter out incorrect point sets and include ones that are correct.

Basic geometry will now come into play as a line is more conveniently represented as a slope and intercept so points at any place on the line can be readily computed. recall that for 2 points, the slope m of a line in a 2 dimensional plane with axes x and y is defined as the change in y over the change in x. The intercept b is computed using one of the points to fit the standard line equation. Fortunately we can use the slope-intercept notation as there are no infinite slopes in time series data where there is a single point per time. Otherwise an polar notation with angle-distance can be used. The distance between the two points is straight-forward from the diagonal of the rectangle formed by the points as given by Pythagorean’s theorem for right triangles.

For 3 points, if we fit the line to 2 points, we can compute the magnitude of the distance from that line to the 3rd point along the y axis. This distance as an error measure is quite arbitrary however because it would be different depending on which 2 points we choose out of 3 possible choices, the 2 additional being: the first and last or the second and third. All of these formula and their examples are given visually:

Closing Price Points Demonstrating Line CalculationsGiven 3 points, a best fit line can be found which is exactly what a linear regression computes. It finds a line which best fits the 3 or more points based on minimizing the error. The error is typically the sum of squared residuals where the residuals are equivalent to the distance measure as given above. The main equation needed is the slope of the line, from which the intercept is computed as per above using the average y and x values.

There are two standard errors given as one is the error of the slope, the other of the intercept. For the slope the error is based on the sum of squared residuals (SSR) computed in part by dividing by the number of points minus 2 where the number of points in this case is 3 making the term disappear so it merely is the square root of sum of squares of the y differences over the sum of squares of x differences. The intercept error is the slope error adjusted by the x values. The slope error is sufficient indicator of the error so this will be applied. It should be noted the SSR, slope and intercept error are assuming that the variance of the points is constant which may not be the case. Assuming the standard slope and intercept errors are normally distributed which also is possibly not true, these would indicate that 68% percent of the points are within plus or minus that error value, where as double the error is 95% and triple the error is 99.7% as per the standard deviation bell curve. So since the variance is not guaranteed nor is a normal distribution, why even use anything derived from the SSR? Since we are not hypothesizing about the distribution but merely comparing the error values to a fixed percentage error, our use case makes the slope error sufficient for selection. The SSR can still be kept as it could also be adapted for the same purpose although with very similar results. However, doing a simple straight average of the expected values verse the actual values is the simplest and still acceptable.

Closing Price Points Demonstrating Linear RegressionIn Python the 3 point variation can be efficiently coded as:

def get_bestfit3(x0, y0, x1, y1, x2, y2):
xbar, ybar = (x0 + x1 + x2) / 3, (y0 + y1 + y2) / 3
xb0, yb0, xb1, yb1, xb2, yb2 =
x0-xbar, y0-ybar, x1-xbar, y1-ybar, x2-xbar, y2-ybar
xs = xb0*xb0+xb1*xb1+xb2*xb2
m = (xb0*yb0+xb1*yb1+xb2*yb2) / xs
b = ybar – m * xbar
ys0, ys1, ys2 =
(y0 – (m * x0 + b)),(y1 – (m * x1 + b)),(y2 – (m * x2 + b))
ys = ys0*ys0+ys1*ys1+ys2*ys2
ser = np.sqrt(ys / xs)
return m, b, ys, ser, ser * np.sqrt((x0*x0+x1*x1+x2*x2)/3)

However, for generality but at the cost of speed, we will implement this for any number of points since a better trend line would be one which could have even more than 3 points given that 3 is merely a minimum:

def get_bestfit(pts):
xbar, ybar = [sum(x) / len (x) for x in zip(*pts)]
def subcalc(x, y):
tx, ty = x – xbar, y – ybar
return tx * ty, tx * tx, x * x
(xy, xs, xx) =
[sum(q) for q in zip(*[subcalc(x, y) for x, y in pts])]
m = xy / xs
b = ybar – m * xbar
ys = sum([np.square(y – (m * x + b)) for x, y in pts])
ser = np.sqrt(ys / ((len(pts) – 2) * xs))
return m, b, ys, ser, ser * np.sqrt(xx / len(pts))

Even more generally, the numpy library has polyfit and poly1d functions which can do exactly this for any polynomial which in this case is a line or degree of 1. We will use it to calculate an average support and average resistance line based on all the local minima and maxima points respectively. The slope and intercept are returned by polyfit and poly1d though poly1d provides easier manipulation and usage and hence is demonstrated where the error returned is merely the sum of squares residual error but for comparison purposes is consistent. Although appropriate for single usage such as the overall average, such a general function might not be as fast as optimized functions written above, hence the exercise of going through writing them. Of course at least 2 points must have been identified to fit this line.

ymin, ymax = [h[x] for x in minimaIdxs], [h[x] for x in maximaIdxs]zmin, zmne, _, _, _ = np.polyfit(minimaIdxs, ymin, 1, full=True) #y=zmin[0]*x+zmin[1]
pmin = np.poly1d(zmin).c
zmax, zmxe, _, _, _ = np.polyfit(maximaIdxs, ymax, 1, full=True) #y=zmax[0]*x+zmax[1]
pmax = np.poly1d(zmax).c
print((pmin, pmax, zmne, zmxe))

If preferring to use the more numerically stable code as per documentation then there is the Polynomial.fit version as well:

p, r = np.polynomial.polynomial.Polynomial.fit
(minimaIdxs, ymin, 1, full=True) #more numerically stable
pmin, zmne = list(reversed(p.convert().coef)), r[0]
p, r = np.polynomial.polynomial.Polynomial.fit
(maximaIdxs, ymax, 1, full=True) #more numerically stable
pmax, zmxe = list(reversed(p.convert().coef)), r[0]

The error cannot be applied consistently for all securities however as absolute values of error are relative to the time period and price range over that time period. So first an appropriate scale should be computed: scale = (hist.Close.max() — hist.Close.min()) / len(hist). (If merely taking square root of the residual error and dividing by (n-2) then dividing by len(hist) here is not necessary.) Then a parameter errpct to the trend lines function will be a simple percentage error such as 0.5%=0.005 where fltpct=scale*errpct. Other nuances should be handled, such as the fact that a slope of 0 is not returned as a coefficient and must be manually filled.

Naive method
This method would simply enumerate through all combinations of 3 points to find the associated error and filter out ones with too great an error value. Of course an O(n³) is not ideal and if the data set were large enough would not even be practically feasible. But the total minima and total maxima points in practice will not be so great that it will not be possible. For example 100 pivot points would be 100*100*100=1,000,000 or a million computations. Obviously to strengthen this with 4 or 5 points would start to become impractical as the order of complexity went to O(n⁴) or O(n⁵). So a different strategy is needed.

def get_trend(Idxs):
trend = []
for x in range(len(Idxs)):
for y in range(x+1, len(Idxs)): for z in range(y+1, len(Idxs)): trend.append(([Idxs[x], Idxs[y], Idxs[z]], get_bestfit3(Idxs[x], h[Idxs[x]], Idxs[y], h[Idxs[y]], Idxs[z], h[Idxs[z]])))
return list(filter(lambda val: val[1][3] mintrend, maxtrend = get_trend(minimaIdxs), get_trend(maximaIdxs)

Sorted Slope method
Fortunately, there are certain properties of the points that can be correlated such as the slope of the line (or the angle from the origin). By doing an O(n²) traversal of all combinations of 2 points, and calculating the slope of the line they form, a list of slopes for each point can be generated. Now sorting a list is worst case complexity ideally of O(n log n) for efficient sorting algorithm like merge sort, and we need to do this over all the n lists for O(n² log n) which will be the complexity of the algorithm. Neighboring points in the sorted slope list can be considered contiguously from 3 to however many points continue to meet the filtering criterion. Once they are matched, the group is removed and the search continues. This part of the algorithm is also O(n²). (The greatest complexity factor is typically the only one to consider to keep the formula simplified so the other 2 O(n²) are not added in.)

This is an approximation algorithm and not exhaustive since sorting by slope is not a guarantee that neighboring slope values will have best fits when there could be great distances between points. However, in practice, it is relatively rare for a scenario to occur where this would happen.

def get_trend_opt(Idxs):
slopes, trend = [], []
for x in range(len(Idxs)): #O(n^2*log n) algorithm
slopes.append([])
for y in range(x+1, len(Idxs)): slope = (h[Idxs[x]] – h[Idxs[y]]) / (Idxs[x] – Idxs[y]) slopes[x].append((slope, y))
for x in range(len(Idxs)):
slopes[x].sort(key=lambda val: val[0])
CurIdxs = [Idxs[x]]
for y in range(0, len(slopes[x])): CurIdxs.append(Idxs[slopes[x][y][1]]) if len(CurIdxs) res = get_bestfit([(p, h[p]) for p in CurIdxs]) if res[3] CurIdxs.sort() if len(CurIdxs) == 3: trend.append((CurIdxs, res)) CurIdxs = list(CurIdxs) else: CurIdxs, trend[-1] = list(CurIdxs), (CurIdxs, res) else: CurIdxs = [CurIdxs[0], CurIdxs[-1]] #restart search
return trend
mintrend, maxtrend =
get_trend_opt(minimaIdxs), get_trend_opt(maximaIdxs)

In fact a lot of matches with 3 points are immediately replaced by 4 or even more points as unsurprisingly trends do occur in real security or index data.

Hough Line Transform method
There are alternative ways to attempt to solve the line finding algorithm, and one comes from image processing where finding lines in images is a common and important task in computer vision. One such method of doing so is a Hough line transform which looks for points across lines on certain angles while scoring them for how many points are found. This algorithm also has limitations. Unfortunately, it uses a lot of memory to keep track of all the histograms and its accuracy is based upon how many angles are tried. The more angles attempted, the slower it becomes however. The memory is based on the diagonal size of the image (with width by height dimensions) and the number of angles (numangles) or ceil(sqrt(width*width+height*height)) * numangles. Scaling of the image is necessary to get good results. In fact the smaller of arctan(2/width) and arctan(2/height) would be the minimum angle to search to find all 3 point possibilities as it is the smallest incremental between a vertical or horizontal line. By fixing the algorithm with a number such as 360*5 possible angles between 90 and -90 degrees, we can use int(np.ceil(2/np.tan(np.pi / (360 * 5)))) to find the maximum image size and merely scale by this amount if the image exceeds the bounds.

We start by adapting our time series data into an image. This requires making the price value discrete which can be done by multiplying by 100 to eliminate the decimal part of the dollar amount. The image only need sized by the length of time and the minimum and maximum price during that time. The whole image will be initialized as black, and white points will be set at the appropriate places where the minima or maxima points occur.

def make_image(Idxs):
max_size = int(np.ceil(2/np.tan(np.pi / (360 * 5)))) #~ m, tested_angles =
hist.Close.min(), np.linspace(-np.pi / 2, np.pi / 2, 360*5)
height = int((hist.Close.max() – m + 0.01) * 100)
mx = min(max_size, height)
scl = 100.0 * mx / height
image = np.zeros((mx, len(hist))) #in rows, columns or y, x
for x in Idxs:
image[int((h[x] – m) * scl), x] = return image, tested_angles, scl, m

The Hough transform is quite straight forward to implement in Python optimized as a point input algorithm. For all angles and all points, the distance to a line perpendicular to the angle which falls through the point is calculated. For each angle, the distance to this perpendicular line is accumulates a point, where different points with identical distances by geometry must lie on a straight line.

Closing Price Points Demonstrating Hough transform accumulation of rho-theta for 2 point line whose origin is based as the first day and minimal price shownThis leads to an O(n*m) algorithm where m is number of angles while memory usage requires m times the diagonal length of the 2-dimensional point space:

def hough_points(pts, width, height, thetas):
diag_len = int(np.ceil(np.sqrt(width * width + height * height)))
rhos = np.linspace(-diag_len, diag_len, diag_len * 2.0)
# Cache some resuable values
cos_t = np.cos(thetas)
sin_t = np.sin(thetas)
num_thetas = len(thetas)
# Hough accumulator array of theta vs rho
accumulator =np.zeros((2 * diag_len, num_thetas), dtype=np.uint64)
# Vote in the hough accumulator
for i in range(len(pts)):
x, y = pts[i]
for t_idx in range(num_thetas): # Calculate rho. diag_len is added for a positive index rho=int(round(x * cos_t[t_idx] + y * sin_t[t_idx])) + diag_len accumulator[rho, t_idx] += 1
return accumulator, thetas, rhos

Hough transform to preserve memory does not return any specific point information, which for visualization purposes can be very useful. So we can calculate the distance for all points to the line, sort it, and take as many points as possible that fit within the error tolerance. The distance of a point to a line is used to determine the correct accumulator to increment. Remember that the slope of the perpendicular slope is the negative reciprocal of the slope so that multiplying them equals -1. A perpendicular line is constructed to the point relative to new point, and the intersection is calculated. The remainder of the proof is derived using the distance between the point and the intersection point. Notice that the numerator is merely the difference between the y and expected y, while the denominator is held constant as we consider only a single line slope. The figure showing the Hough transform also gives the explicit formula.

def find_line_pts(Idxs, x0, y0, x1, y1):
s = (y0 – y1) / (x0 – x1)
i, dnm = y0 – s * x0, np.sqrt(1 + s*s)
dist = [(np.abs(i+s*x-h[x])/dnm, x) for x in Idxs]
dist.sort(key=lambda val: val[0])
pts, res = [], None
for x in range(len(dist)):
pts.append((dist[x][1], h[dist[x][1]]))
if len(pts) r = get_bestfit(pts)
if r[3] > fltpct: pts = pts[:-1] break
res = r
pts = [x for x, _ in pts]
pts.sort()
return pts, res

We can also use the Hough line transform function of the scikit-image library called hough_line. It should be noted that the OpenCV (Computer Vision) library which is also available in Python, also has the same functions. This could also be hand implemented in a short amount of code since the algorithm is not particularly complex but the libraries are optimized and faster. Notice that only 2 points in the accumulator is enough for the algorithm to return them so we would filter for 3 or more associated points. Useful details are that it reverses the axes and has a threshold parameter which is greater but not equal to. There are some small detailed differences in how it functions internally especially with regard to the selection of local maxima in the accumulator so the results will not be exactly identical. This leads to two different methods for performing the calculation (one point-optimized, the other transforming points to image for the library):

def houghpt(Idxs):
max_size = int(np.ceil(2/np.tan(np.pi / (360 * 5)))) #~ m, tested_angles =
hist.Close.min(), np.linspace(-np.pi / 2, np.pi / 2, 360*5)
height = int((hist.Close.max() – m + 1) * 100)
mx = min(max_size, height)
scl = 100.0 * mx / height
acc, theta, d = hough_points(
[(x, int((h[x] – m) * scl)) for x in Idxs], mx, len(hist),
np.linspace(-np.pi / 2, np.pi / 2, 360*5))
origin, lines = np.array((0, len(hist))), []
for x, y in np.argwhere(acc >= 3):
dist, angle = d[x], theta[y]
y0, y1 = (dist – origin * np.cos(angle)) / np.sin(angle)
y0, y1 = y0 / scl + m, y1 / scl + m
pts, res = find_line_pts(Idxs, 0, y0, len(hist), y1)
if len(pts) >= 3: lines.append((pts, res))
return lines
mintrend, maxtrend = houghpt(minimaIdxs), houghpt(maximaIdxs)def hough(Idxs): #pip install scikit-image
image, tested_angles, scl, m = make_image(Idxs)
from skimage.transform import hough_line, hough_line_peaks
h, theta, d = hough_line(image, theta=tested_angles)
origin, lines = np.array((0, image.shape[1])), []
for pts, angle, dist in zip(*hough_line_peaks(h, theta, d, threshold=2)):
y0, y1 = (dist – origin * np.cos(angle)) / np.sin(angle)
y0, y1 = y0 / scl + m, y1 / scl + m
pts, res = find_line_pts(Idxs, 0, y0, image.shape[1], y1)
if len(pts) >= 3: lines.append((pts, res))
return lines
mintrend, maxtrend = hough(minimaIdxs), hough(maximaIdxs)

Probabilistic Hough Line Transform method
The precise angles used and high processing and memory requirements required for a normal Hough line identification makes it somewhat impractical for this task with a small number of points. The greater the number of points, the better it could work. However this prior method was more of an exercise in precision. In practice a probabilistic Hough line transform is used where the search is done in a random manner, and parameters are used to filter the results including a threshold for the number of points. Another library function probabilistic_hough_line works for this purpose:

def prob_hough(Idxs): #pip install scikit-image
image, tested_angles, scl, m = make_image(Idxs)
from skimage.transform import probabilistic_hough_line
lines = []
for x in range(hough_prob_iter):
lines.append(probabilistic_hough_line(image, threshold=2, theta=tested_angles, line_length=0, line_gap=int(np.ceil(np.sqrt( np.square(image.shape[0]) + np.square(image.shape[1]))))))
l = []
for (x0, y0), (x1, y1) in lines:
if x0 == x1: continue
if x1 y0, y1 = y0 / scl + m, y1 / scl + m
pts, res = find_line_pts(Idxs, x0, y0, x1, y1)
if len(pts) >= 3: l.append((pts, res))
return l
mintrend, maxtrend = prob_hough(minimaIdxs), prob_hough(maximaIdxs)

There are still only 2 points of the line returned here, but the threshold is built in and done for us. Unfortunately, the results will be different every time it is run due to the nature of a probabilistic algorithm having randomness. This in general makes this method not particularly well suited for finding all trend lines. Since it can be run multiple times, this could be done until a certain number of runs has occurred or a certain number of lines is found to greatly increase the probability of finding all of the lines. Thus, the parameter hough_prob_iter specifies the number of iterations to run it such as 10 iterations to increase likelihood of a sufficient number of lines.

Note that better than hard-coding 100 for a scale of pennies or 1/0.01 since dollars are in use, it would be better to have a hough_scale parameter.

We now have 5 different methods for trying to find trend lines.

Unfortunately, just because local minima or maxima points are on a line, they may not always yield the best trend line. This is because what is considered best is often subjective. However there are ways this can be turned into an objective definition. One such way is to see how often the price crosses over the trend line, or the area made by the trend line and the price data across its earliest and latest point. The error percentage is only used to determine what may or may not be a trend line, but not how useful that trend line may or may not be. With enough data, lines can be found years apart which are unlikely to be useful whether or not they are coincidental. Typically trend lines are drawn which are containing the signal on one side or another, so the area based approach seems at least reasonable. A Riemann sum provides an integration technique which can be used to calculate this area. In fact, the application is trivial as its merely subtracting the two and summing all values less or greater than 0. Specifically this will be a right Riemann sum which uses the function value at each time point as opposed to a mid-point or the preceding, maximal or minimal one. Finally dividing it by the number of days gives a measure of the wrong side of the trend per day and the best trend lines can be selected.

Closing Price with Resistance and Areadef measure_area(trendline, isMin):
base = trendline[0][0]
m, b, ser =
trendline[1][0], trendline[1][1], h[base:trendline[0][-1]+1]
return sum([max(0, (m * (x+base) + b) – y if isMin else y – (m * (x+base) + b)) for x, y in enumerate(ser)]) / len(ser)
mintrend = [(pts, (res[0], res[1], res[2], res[3], res[4], measure_area((pts, res), True))) for pts, res in mintrend]
maxtrend = [(pts, (res[0], res[1], res[2], res[3], res[4], measure_area((pts, res), False))) for pts, res in maxtrend]
mintrend.sort(key=lambda val: val[1][5])
maxtrend.sort(key=lambda val: val[1][5])
print((mintrend[:5], maxtrend[:5]))

Another issue, is that these algorithms can be slow even the best of them when efficiently implemented for a variety of reasons. Further optimization is often possible especially by eliminating pandas and moving towards lists and numpy arrays which have fast C implementations.

The real question here though, is how exactly are the trend lines being used? Are they meant to be short term or long term, current or historical? This is important enough that the idea of calculating trend lines over windows must be dealt with. How exactly can a window be made here?

The first idea is to take a window of say a single quarter or a year as too short of a window will not yield enough pivot points to be useful. Then the window must be rolled through all the data, however, not in such a way that its inefficient or duplicating the work, nor is it missing possible trend lines. Once the window size is decided upon, the work can begin by searching with double the window size. Doubling the window size ensures that trend lines across the borders of the windows will be found. Then the double window size is searched across all of the data in window sized increments. So if your window is a year, you search the past two years, then search from 3 years prior to the past year, then from 4 years prior to 2 years prior, etc. One part of the window will end up getting searched twice and many identical or extensions of identical trends will appear. This will be much more visually pleasing as well, since it would be too disruptive to the plot to have trend lines drawn too distantly.

So finally, the trend lines found must be intelligently merged back together to get the final set. It can be done by going through all trend lines containing each given point and applying the previously discussed slope sorting method. This is necessary for the Hough transform methods anyway as they can yield redundant point sets.

def merge_lines(Idxs, trend):
for x in Idxs:
l = []
for i, (p, r) in enumerate(trend): if x in p: l.append((r[0], i))
l.sort(key=lambda val: val[0])
if len(l) > 1: CurIdxs = list(trend[l[0][1]][0])
for (s, i) in l[1:]: CurIdxs += trend[i][0] CurIdxs = list(dict.fromkeys(CurIdxs)) CurIdxs.sort() res = get_bestfit([(p, h[p]) for p in CurIdxs]) if res[3] ([], None), (CurIdxs, res), list(CurIdxs) else: CurIdxs = list(trend[i][0]) #restart search from here
return list(filter(lambda val: val[0] != [], trend))
mintrend, maxtrend = merge_lines(minimaIdxs, mintrend), merge_lines(maximaIdxs, maxtrend)

The main take away from this discussion is to objective your idea of best, knowing that not all people will necessarily agree on any one objective definition. It helps to not take any shortcuts during the brainstorming stage to create a comprehensive and proper method.

The merging process done above could in fact find trend lines which happen to have points in windows very far apart which would extensively lengthen them, something which may or may not be desired. For display purposes however, re-windowing the trend lines found can be ideal. Regardless of how they are merged they can be split back apart by their windows. Generally we want to draw the trend lines only into the future for a limited range to avoid cluttering the plot.

Closing Price with best 2 Support/Resistance Trend Lines by errorNotice in the figure that the bottom support line looks from a financial perspective as correct and useful as there are no occasions where the price dropped below it. The other support and resistance lines chose only semi-prominent peaks and are legitimate per definition of trend line, but hardly what would be considered useful. So the best trend lines with lowest error, are not necessarily the best support and resistance lines. Again the idea of measuring area formed by the above the line to stock price for resistance, and below the line to the stock price for support seems to be the key.

Closing Price with best 2 Support/Resistance Trend Lines by smallest wrong side area per dayThese trend lines are quite good for the ranges of points they contain. Non-prominent local extrema are still going to sometimes yield trend lines of limited use like the flat resistance line. Also very short term trends like the very steep support line might not be relevant so far into the future, and further heuristics could be used to try to eliminate them though it could be considered useful as it was previously in that short term. However after its third point, as that was a recovery period, the market reversed course and by the time it was identified its usefulness was already outdated. The other best support and resistance line are absolutely what an analyst would hope to see. Simply not drawing out into the future might also help but for the most recent window, the whole purpose is to project towards the future.

Lastly, instead of showing the 2 best support/resistance lines for the whole period, the final result can be re-windowed for display and extrapolation purposes and the best 2 from each window can be used.

The trend lines can be used to extrapolate feature data for machine learning such as using an LSTM (Long-Short Term Memory) or GRU (Gated Recurrent Unit) neural network even as part of a GAN (Generative Adversarial Network) for predicting the time series, to mention some realistic modern ideas. Perhaps they will first be reduced with the other features by use of PCA (Principal Component Analysis) or auto-encoders.

However, to prevent data leak, any trend line found can only be used to predict points in the future after the last point on the trend line. Any time at or before the final point, would be a classic example of data leak where future values are used to leak values into the past. Even selecting best trend lines must be carefully done only using data prior to the final point, although visualization does not have such a limitation. Caution must always be careful here as even typical programming bugs such as being off by 1 in computations are fatal mistakes when dealing with a precision task like time series prediction. Remember the detail of the numerical differentiation which does look ahead as well when dealing with the first day or two depending on accuracy, after a trend has occurred. It could be better to assume the trend starts a day or two after the last minima or maxima point to avoid this kind of leak depending on how the trend is being used.

There is a GitHub repository for trendln containing all of the code for perusal. This is also available as a PyPI package trendln which can be installed and used readily.

Dependencies include numpy. For plotting matplotlib and pandas are required. For image-based Hough trend lines and Probabilistic Hough trend lines, scikit-image is required. For numerical differentiation, findiff is required. The code which generated all the images in this article is also included as reference material only.

We have analyzed and demonstrated how to use different methods including numerical differentiation to find pivot points. Additionally we have shown how to measure error in lines as well as use various trend line identification methods including using a sorted slope, Hough line transforms and the probabilistic Hough line transform. Now correct trend lines for support and resistance can be collected and studied. They can be further applied for more complex analysis or used for forecasting movement in security pricing.

Furthermore, there was an incentive to share this given that it combines a lot of information into one resource with a thorough analysis of each aspect of the task while most claimed solutions have not even met the basic definition of a trend line.

If you enjoyed this article kindly let me know and feel free to ask questions. Stay tuned for shape identification techniques as was alluded to.