Curve Fitting is a Gram-Schmidt Reduction
One of the things you can do with time series data is to fit a curve (say trend, or average, or natural cycles) and then analyze the residue. The subject is studied to death out there, but I discovered a neat trick that I would like to share.
Assume we have a function \(f(x)\), but we want to express it in terms of a (finite) collection of other functions \(f_1(x),\ldots,f_n(x)\). In the discrete version of the same problem we have a set of data points \((x_1,\ldots,x_N)\) and we want a collection of data points \((y^{(i)}_1,\ldots,y^{(i)}_{M_i})\) to fit. If without loss of generality, we assume each \(M_i\) is the same as \(N\), we can use the Gram-Schmidt Process to get the data fit. The idea is simple enough to implement, and there is already a nice paper on the subject.
Code
The only non-trivial part of the code is the Gram-Schmidt reduction:
(defun gram-schmidt (xs &optional tail)
(labels ((dot (a b) (reduce #'+ (mapcar #'* a b))))
(if (null (cdr xs))
(cons (car xs) tail)
(destructuring-bind
(y . ys) xs
(gram-schmidt
(mapcar (lambda (z)
(let ((c (/ (dot y z) (dot y y))))
(mapcar (lambda (i j) (- i (* c j))) z y)))
ys)
(cons y tail))))))
GRAM-SCHMIDTThe data I am going to use is from Federal Reserve. I am going to look at the total production index from March 1986 to February 2015. The data need to be extracted from the file I downloaded. I used
grep grep B50001 ip_sa.txt | grep -v : | awk '{ for(i=3;i<15;i++) { print $i; } }' > total_index.csvto extract the pieces I needed. I am going to load the data into a
variable called *raw-data*
(defvar *raw-data*
(with-open-file (in "total_index.csv" :direction :input)
(do* ((item (read-line in nil) (read-line in nil))
(res (list item) (cons item res)))
((null item) (mapcar #'read-from-string (reverse (cdr res)))))))
*RAW-DATA*I am going to fit a degree 5 polynomial with the degree 4 monomial missing onto the data.
(defun residues (raw degrees)
(labels ((poly (k n) (loop for i from 0 below n collect (expt i k))))
(let ((n (length raw)))
(car (gram-schmidt (reverse (cons raw (mapcar (lambda (x) (poly x n)) degrees))))))))
RESIDUES
(defvar *residue* (residues *raw-data* '(0 1 2 3 5)))
*RESIDUE*Let us see the results:
