``````# Copyright 2007, 2009, 2010, 2011 Kevin Ryde

# This file is part of Chart.
#
# Chart is free software; you can redistribute it and/or modify it under the
# Foundation; either version 3, or (at your option) any later version.
#
# Chart is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along
# with Chart.  If not, see <http://www.gnu.org/licenses/>.

package App::Chart::Series::Derived::LinRegStderr;
use 5.010;
use strict;
use warnings;
use Carp;
use List::Util qw(min max);
use Locale::TextDomain ('App-Chart');

use base 'App::Chart::Series::Indicator';
use App::Chart::Series::Derived::SMA;

# http://mathworld.wolfram.com/LeastSquaresFitting.html
#     Showing how to calculate the variance in the error amounts (e[i])
#     without calculating a,b parameters.
#
#
#
#               Sum (y - (a+b*x))^2
# stderr = sqrt -------------------
#                       N^2
#
#               Sum ((y - My) - b*x)^2
#        = sqrt ----------------------     where My = Mean(Y)
#                       N^2
#
#               Sum (y - My)^2 - 2*My*b*x + b^2*x^2
#        = sqrt -----------------------------------
#                               N^2
#
#                   (y - My)^2       2*My*b*x - b^2*x^2
#        = sqrt Sum ---------- - Sum ------------------
#                    N^2                    N^2
#
# The first term is Variance(Y), and with b = Mean(X*Y) - Mean(X) * Mean(Y)
#
#                            2*M(y)*Sum(x) - M(X*Y)*Sum(x^2)/M(X)*M(Y)
#        = sqrt Var(Y) - b * -----------------------------------------
#                                               N^2
#
# ....

sub longname  { __('Linear Regression Stderr') }
sub shortname { __('Linreg Stderr') }
sub manual    { __p('manual-node','Linear Regression Standard Error') }

use constant
{ type       => 'indicator',
units      => 'price',
priority   => -10,
minimum    => 0,
parameter_info => [ { name    => __('Days'),
key     => 'linregstderr_days',
type    => 'integer',
minimum => 1,
default => 20 } ],
};

sub new {
my (\$class, \$parent, \$N) = @_;

\$N //= parameter_info()->->{'default'};
(\$N > 0) || croak "LinRegStderr bad N: \$N";

return \$class->SUPER::new
(parent     => \$parent,
N          => \$N,
parameters => [ \$N ],
arrays     => { values => [] },
array_aliases => { });
}
*warmup_count = \&App::Chart::Series::Derived::SMA::warmup_count; # \$N-1

# Return the factor for 1/Variance(X) arising in the stderr formula.  With
# X values -1.5,-0.5,0.5, 1.5 this means
#
#                         1
#      --------------------------------------------
#      n * ((-1.5)^2 + (-0.5)^2 + (0.5)^2 + (1.5)^2)
#
# and the denominator here is (n^4 - n^2)/12.  See devel/ema-omitted.pl for
# checking this.  Same func in RSquared.pm too.
#
sub xfactor {
my (\$N) = @_;
if (\$N < 2) { return 1; }
return 12.0 / (\$N*\$N * (\$N*\$N - 1));
}

sub proc {
my (\$class, \$N) = @_;

# @array,\$pos is a  circular list of the last \$N many Y values.
#
# \$count is the number of values in @array.
#
# \$y_total is the sum of the values in @array.
#
# \$yy_total is the sum of the squares of the values in @array.
#
# \$x_total is the sum of the X values corresponding to the values in
# @array, which means 1+2+...+\$count, which is a triangular number.  This
# varies until \$count==\$N is reached, and is then a constant.
#
# \$x_factor is the variance factor in the denominator, see `xfactor' above.
# This varies until \$count==\$N is reached, and is then a constant.
#
# \$xy_total is the sum of the product x*y for each x,y pair, which means
# 1*y + 2*y + ... + \$count*y[\$count].  When shifting out an old Y,
# the weighting of each y in the sum is reduced by subtracting \$y_total.

my @array;
my \$pos = \$N - 1;  # initial extends
my \$y_total = 0;
my \$yy_total = 0;
my \$x_factor = 0;
my \$xy_total = 0;

my \$count = 0;
return sub {
my (\$y) = @_;

if (\$count >= \$N) {
# drop oldest point
my \$prev = \$array[\$pos];
\$yy_total -= \$prev ** 2;
\$xy_total += (\$count-1)*0.5 * \$prev;
\$y_total -= \$prev;
} else {
# gaining a point
\$count++;
\$x_factor = xfactor(\$count);
\$xy_total += \$y_total * 0.5; # adj so 0.5 less each x
}
\$array[\$pos] = \$y;
if (++\$pos >= \$N) { \$pos = 0; }

# shift xy products onto 1 less x each
\$xy_total -= \$y_total;

\$y_total += \$y;
\$yy_total += \$y*\$y;
\$xy_total += \$y * (\$count-1)*0.5;

return (sqrt(max (0,
\$count*\$yy_total - \$y_total*\$y_total  # Var(Y)
- ((\$count * \$xy_total)**2   # (Covar(X,Y))^2
* \$x_factor)))       # 1 / (X variance)
/ \$count);
};
}

1;
__END__

#
# App::Chart::Series::Derived::LinRegStderr -- linear regression standard error
#