| 1 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 2 | // FILE: interpolate.tol |
|---|
| 3 | // PURPOSE: Interpolation functions |
|---|
| 4 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 5 | |
|---|
| 6 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 7 | Matrix FillMissingInterpMat(Matrix X, Matrix Y, Text method_, Real order) |
|---|
| 8 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 9 | { |
|---|
| 10 | Text method = If(method_=="","AkimaSpline",method_); |
|---|
| 11 | Set idx = Range(1,Rows(Y),1); |
|---|
| 12 | Set idxKnown = Select(idx, Real(Real r) |
|---|
| 13 | { |
|---|
| 14 | !IsUnknown(MatDat(Y,r,1)) |
|---|
| 15 | }); |
|---|
| 16 | Matrix y = SubRow(Y,idxKnown); |
|---|
| 17 | Matrix x = SubRow(X,idxKnown); |
|---|
| 18 | If(method_==ToLower(method_), |
|---|
| 19 | { |
|---|
| 20 | //GSL interpolation |
|---|
| 21 | Code interpolation = gsl_interp(method,x,y); |
|---|
| 22 | SetCol(EvalSet(idx,Real(Real r) |
|---|
| 23 | { |
|---|
| 24 | Real x_ = MatDat(X,r,1); |
|---|
| 25 | interpolation(0,x_) |
|---|
| 26 | })) |
|---|
| 27 | }, |
|---|
| 28 | { |
|---|
| 29 | //AlgLib interpolation |
|---|
| 30 | Code interpolation = AlgLib.Interp.Scalar(method,x,y); |
|---|
| 31 | SetCol(EvalSet(idx,Real(Real r) |
|---|
| 32 | { |
|---|
| 33 | Real x_ = MatDat(X,r,1); |
|---|
| 34 | interpolation(x_) |
|---|
| 35 | })) |
|---|
| 36 | }) |
|---|
| 37 | }; |
|---|
| 38 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 39 | PutDescription("Fill missing values of Y by interpolation of function Y(X)." |
|---|
| 40 | "The used interpolation engine will be gsl_interp if method is one of " |
|---|
| 41 | "these:\n" |
|---|
| 42 | " linear, polynomial, cspline, cspline_periodic, akima,akima_periodic.\n" |
|---|
| 43 | "The used interpolation engine will be AlgLib.Interp.Scalar if method is " |
|---|
| 44 | "one of these:\n" |
|---|
| 45 | " BarycentricRational, LinearSpline, CubicSpline, AkimaSpline, " |
|---|
| 46 | "SplineLeastSquares, ChebyshevLeastSquares, PolynomialLeastSquares\n" |
|---|
| 47 | "Default method is AkimaSpline", |
|---|
| 48 | FillMissingInterpMat); |
|---|
| 49 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 50 | |
|---|
| 51 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 52 | Serie FillAllMissingInterp(Serie y, Text method_, Real order) |
|---|
| 53 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 54 | { |
|---|
| 55 | Text method = If(method_=="","AkimaSpline",method_); |
|---|
| 56 | Matrix Y = Tra(SerMat(y)); |
|---|
| 57 | Matrix X = SetCol(Range(1,Rows(Y),1)); |
|---|
| 58 | Matrix Y_ = FillMissingInterpMat(X, Y, method, order); |
|---|
| 59 | MatSerSet(Tra(Y_),Dating(y),First(y))[1] |
|---|
| 60 | }; |
|---|
| 61 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 62 | PutDescription("Fill all missing values of a serie by interpolation of function " |
|---|
| 63 | "y(t) of natural index in discrete time. Missing intervals with length " |
|---|
| 64 | "greater than order will not be filled.\n" |
|---|
| 65 | "Interpolation method is the same used in FillMissingInterpMat", |
|---|
| 66 | FillAllMissingInterp); |
|---|
| 67 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 68 | |
|---|
| 69 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 70 | Serie FillMissingInterp(Serie y, Text method_, Real order) |
|---|
| 71 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 72 | { |
|---|
| 73 | Text method = If(method_=="","AkimaSpline",method_); |
|---|
| 74 | Serie y1 = !IsUnknown(y); |
|---|
| 75 | Serie yo = Or(1,y1) << (Expand((1-B^order)/(1-B),order):y1); |
|---|
| 76 | TimeSet CtYo = SerTms(yo); |
|---|
| 77 | Serie yok = DatCh(y, CtYo, FirstS); |
|---|
| 78 | Matrix Y = Tra(SerMat(yok)); |
|---|
| 79 | Matrix X = SetCol(Range(1,Rows(Y),1)); |
|---|
| 80 | Matrix Y_ = FillMissingInterpMat(X, Y, method, order); |
|---|
| 81 | Serie yok_ = MatSerSet(Tra(Y_),CtYo,First(yok))[1]; |
|---|
| 82 | InvCh(yok_, y) |
|---|
| 83 | }; |
|---|
| 84 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 85 | PutDescription("Fill missing values of a serie by interpolation of function " |
|---|
| 86 | "y(t) of natural index in discrete time. Missing intervals with length " |
|---|
| 87 | "greater than order will not be filled.\n" |
|---|
| 88 | "Interpolation method is the same used in FillMissingInterpMat", |
|---|
| 89 | FillMissingInterp); |
|---|
| 90 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 91 | |
|---|
| 92 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 93 | Serie InvChInterp(Serie S1, TimeSet dating, Text method_, Real order) |
|---|
| 94 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 95 | { |
|---|
| 96 | Serie serExt = InvCh(S1, CalInd(W,dating)*?); |
|---|
| 97 | Serie aux = FillAllMissingInterp(serExt, method_, order); |
|---|
| 98 | Eval(ToName(Name(S1))+".Interp=aux") |
|---|
| 99 | }; |
|---|
| 100 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 101 | PutDescription("Extends time series to a more fine dating by interpolation.\n" |
|---|
| 102 | "Interpolation method is AlgLib.Interp.Scalar. Default method is AkimaSpline", |
|---|
| 103 | InvChInterp); |
|---|
| 104 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 105 | |
|---|
| 106 | |
|---|
