Statistics

petros ypsilantis
May 2022

1 Introduction
Sum of Independent RVs
It is always true (even when X1 , X2 , ....) that

E(X1 + X2 + ... + Xn ) = E(X1 ) + E(X2 ) + ... + E(Xn ) (1)

if X1 , X2 , .... independent we have E(X1 X2 ...Xn ) = E(X1 )E(X2 )...E(Xn )

and V ar(X1 + X2 + ... + Xn ) = V ar(X1 ) + V ar(X2 ) + ... + V ar(Xn )
Also, we have,

V ar(a1 X1 +a2 X2 +...+an Xn ) = a21 V ar(X1 )+a22 V ar(X2 )+...+a2n V ar(Xn ) (2)

if X1 , X2 , .... not independent then we have

n
X XX
V ar(X1 + X2 + ... + Xn ) = V ar(Xi ) + 2 Cov(Xi , Xj ) (3)
i=1 i<j
Xn XX q
= V ar(Xi ) + 2 ρi,j V ar(Xi )V ar(Xj )
i=1 i<j
(4)

For example, if we have 2 stocks A and B, and the variance of stock’s A

2 2
return is σA , for B is σB and they also have correlation ρ, then we have that

V ar(rA − hrB ) = V ar(rA ) + h2 V ar(rB ) − 2hρσA σB (5)

Laws of Large Numbers:
Let X1 , X2 , ... is a sequence of IID RVs, each having finite mean µ and
variance σ 2 . Let Sn := X1 + X2 + ... + Xn so that E(Sn ) = nµ, V ar(Sn ) = nσ 2 .
The weak law of large numbers is that for ϵ > 0,

P (|Sn /n − µ| > 0) → 0, as n → ∞ (6)

The strong law of large numbers is that

1
 P (Sn /n → µ) = 1 (7)
Central Limit Theorem
Let’s assume that X1 , ... as above. The CLT states that the standardized
version Sn∗ of Sn as
Sn − nµ
Sn∗ = √ (8)
σ n
so that E(Sn∗ ) = 0 and V ar(Sn∗ ) = 1. Then for large n, Sn∗ is approximately
standard normal. Another way to see this if X̄n is the sample mean √
of the first n
samples, then the limiting form of the distribution Z = limn→∞ n(X̄σn −µ) , is a
standard normal distribution. For example, suppose that a sample is obtained
containing many observations, each observation being randomly generated in
a way that does not depend on the values of the other observations, and that
the arithmetic mean of the observed values is computed. If this procedure
is performed many times, the central limit theorem says that the probability
distribution of the average will closely approximate a normal distribution. A
simple example of this is that if one flips a coin many times, the probability of
getting a given number of heads will approach a normal distribution, with the
mean equal to half the total number of flips. At the limit of an infinite number
of flips, it will equal a normal distribution.
The central limit theorem has several variants. In its common form, the
random variables must be identically distributed. In variants, convergence of
the mean to the normal distribution also occurs for non-identical distributions
or for non-independent observations, if they comply with certain conditions.
Maximum Likelihood Estimation
Suppose a RV X has pdf f (x|θ). Given the observed value x of X, the likeli-
hood of θ is defined as lik(θ) = f (x|θ). In the case of multiple observations, i.e.,
when x = (x1 , x2 , ..., xn ) is a vector of observed values of X1 , ..., Xn we assume,
unless otherwise stated, that X1 , ..., Xn are IIDs; in this case f (x1 , ..., xn |θ) is
the product of marginals.
n
Y
lik(θ) = f (x1 , ..., xn |θ) = f (xi |θ) (9)
i=1

It makes intuitive sense to estimate θ by whatever value gives the greatest

likelihood to the observed data. Thus the maximum likelihood estimate θ̂(x) of
θ is defined as the value of θ that maximizes the likelihood. Then θ̂(x) is called
the maximum likelihood estimator (MLE) of θ.
Of course, the maximum likelihood estimator need not exist, but in many
examples it does. In practice, we usually find the MLE by maximizing logf (x|θ),
which is known as the loglikelihood
MLE Properties

1. If it exists it is a function of sufficient statistic. The informal notion

of a sufficient statistic T = T (X1 , X2 , ...Xn ) is that it summarizes all

2
 information in {X1 , ..., Xn }. Formally, the statistic T = T (X) is said to
be sufficient for θ if, for each t, Pθ (X ∈ .|T (X) = t) does not depend on
θ. I.e, the conditional distribution of X1 , ..., Xn given T (X) = t doesn’t
involve θ. Thus to know more about x that T (x) = t is of no additional
help in making any inference about θ.
2. MLE is asymptotically Unbiased Estimator: E(θ̂) → θ as n → ∞. Now
we need to understand how fast it converges to the parameter.
3. Rao Blackwell Theorem (Mean square Error): A good estimator should
take values close to the true value of the parameter it is attempting to
estimate. If θ̂ is an unbiased estimator of θ then E(θ̂ − θ)2 is the variance
of θ̂. If θ̂ is a biased estimator then E(θ̂ − θ)2 is no longer the variance of
θ̂ but it is still a useful measure of the mean square error (MSE) of θ̂. The
theorem says that if we want an estimator with small MSE we can confine
our search to estimators which are functions of the sufficient statistic. Rao
Blackwell Theorem: Let θ̂ be an estimator of θ with E(θˆ2 ) < ∞ for all θ.
Suppose that T is sufficient for θ, and let θ∗ = E(θ̂|T ). Then, for all θ,

E(θ∗ − θ)2 <= E(θ̂ − θ)2 (10)

The inequality is strict unless θ̂ is a function of T .
4. Consistency We say that θ̂ is consistent if

P (|θ̂ − θ| > ϵ) → 0as n → ∞ (11)

5. Asymptotically Efficient It can be be shown that V ar(θ̂) >= 1/nI(θ)

for any unbiased estimate θ̂, where 1/nI(θ) is called the Cramer-Rao lower
bound. We say that θ̂ is asymptotically efficient if

var(θ̂)
limn→∞ =1 (12)
[1/nI(θ)]
6. Under certain conditions it can be shown √ that the MLE converges in dis-
tribution to a normal distribution. It is n-consistent and asymptotically
efficient, meaning that it reaches the Cramer Rao Lower bound.
√
n(θ̂ − θ) → N (0, I −1 ), (13)

where I is the fisher information matrix.

∂ 2 ln fθ0 (Xt )
Ijk = E[ ] (14)
∂θj ∂θk
In Particular it means√that the bias of the MLE is equal to
zero up to the order 1/ n. The fact that MLE asymptotically follows

3
 Normal distribution can be used to estimate confidence intervals for the
parameter. If we have a single parameter then we have,

1.96
θ̂ + − q (15)
nI(θ̂)

1
Pn ∂ 2 ln fθ(Xi )
where I(θ̂) = n i=1 ∂2θ

2 The Lasso
The lasso estimate is defined by
N
X p
X
β̂ lasso = argminβ (yi − β0 − xij βj )2 (16)
i=1 j=1
Pp
subject to j=1 |βj | <= t. We can write the lasso problem in the equivalent
Lagrangian form,
N
X p
X p
X
β̂ lasso = argminβ { (yi − β0 − xij βj )2 + λ |βj |} (17)
i=1 j=1 j=1

The constrain makes the solution non-linear in the yi , and there is no closed
form expression as in the ridge regression. Computing the lasso solution is
a quadratic programming problem, although there are efficient algorithms for
solving the problem.

2.1 Lasso Interpretation

The below figure depicts the lasso (left) and ridge regression (right) when there
are only 2 parameters. The residual sum of squares has elliptical contours,
centered at the full least squared estimate. The constraint region for ridge
regression is the disk (β12 = β22 ) <= t, while that for lasso is the diamond
|β1 | + |β2 | <= t. Both methods find the first point where the elliptical contours
hit the constraint region. Unlike the disk, the diamond has corners; if the
solution occurs at a corner, then it has one parameter j equal to zero. When
p > 0, the diamond becomes a rhomboid, and has many¯ corners, flat edges and
faces; their are many opportunities for the estimated parameter to be zero.

4
 Another reasoning on how L1-regularization is more likely to create 0-weights.
Consider a model consisting of the weights (w1 , w2 , ..., wm ). P
With L1 regular-
ization, you penalize the model by a loss function L1 (w) = Pi wi . With L2-
regularization, you penalize the model by a loss L2 (w) = 21 i wi2 . If we use
gradient descent, you will iteratively make the weights change in the opposite
direction of the gradient with a step size η multiplied with the gradient. This
mean that a more steep gradient will make us take a larger step, while a more
flat gradient will make us take a smaller step. Let us look at the gradients
(subgradient in case of L2): dLdw
1 (w)
= sign(w), where sign(w) = ( |w w1
, w1 , ...)
1 | |w1 |

and dLdw
2 (w)
=w
Notice that for L1 , the gradient is either 1 or −1, except for when w1 = 0.
That means that L1-regularization will move any weight towards 0 with the
same step size, regardless the weight’s value. In contrast, you can see that
the L2 gradient is linearly decreasing towards 0 as the weight goes towards 0.
Therefore, L2-regularization will also move any weight towards 0, but it will
take smaller and smaller steps as a weight approaches 0.

5
 Try to think you start with a model with w1 = 5 and using η = 12 . In the
following picture you can see how gradient descent using gradient descent using
L1 regularization makes 10 of the updates. w1 := w1 − η dLdw 1 (w)
= w1 − 12 · 1,
until reaching a model with w1 = 0

In contrast, with L2-regularization where η = 12 , the gradient w1 , causing

every step to be only halfway towards 0. That is, we make the update w1 :=

6
w1 − η · dLdw
2 (w)
= w1 − 12 · w1 . Therefore, the model never reaches a weight of 0,
regardless of how many steps we take:

3 Regularization in Logistic Regression

Does the regularization in LR always results in better fit and better generaliza-
tion?
Answer: Regularization does not improve the performance on the data set
that the algorithm used to learn the model parameters. However it can improve
the generalization performance.

4 Bayesian Linear Regression

The frequentist, or classical approach to multiple linear regression assumes a
model of the form
p
X
f (X) = β0 + Xj βj + ϵ = β T X + ϵ, (18)
j=1

where ϵ ∼ N (0, σ 2 ).
The syntax for a linear regression in a Bayesian framework look like this:

y ∼ N (β T X, σ 2 I) (19)
What is the benefit of doing this?

1. Prior Distribution: If we have any prior knowledge about the parameters

β, then we can choose prior distributions that reflect this. If we don’t we
can choose non-informative priors.

7
 2. Posterior Distributions: In the frequentist approach, β̂ was only a single
point estimate. In the BAyesian formulation we receive an entire prob-
ability distribution that characterises our uncertainty on the different β
coefficients. The immediate benefit of this is that after taking account any
data we can quantify our uncertainty in the β parameters via the variance
of this posterior distribution.

while the above formula for the Bayesian approach may appear succinct, it
doesn’t really give us much clue as to how specify a model and sample from it
using MCMC.
Let’s start with Generalized Linear Models?
Here we discuss the concept of Generalized Linear Models as we’ll be using
these to formulate our model in PyMC3. A GLM is a flexible mechanism for
extending ordinary LR to more general form of regression, including logistic re-
gression and poisson regression (used for count data) as well as linear regression
itself.
GLMs allow for response variables to have error distributions other than
Normal. The linear model is related to the response/outcome, y via a link
function and is assumed to be generated from a statistical distribution from the
exponential distribution family. This family of distributions encompasses many
common distributions including normal, gamma, beta, chi-squared, Bernoulli,
Poisson and others.
The mean of this distribution µ depends on X via the following relation:

E(y) = µ = g −1 (Xβ), (20)

where g is the link function. The variance in often some function V , of the
mean,

V ar(y) = V (E(y)) = V (g −1 (Xβ)) (21)

In the frequentist setting, as with ordinary linear regression above, the unknown
coefficients are estimated via a maximum likelihood approach. I’m not going
to discuss GLMs in depth here as they are not the focus of the article. We
are interested in them because we will be using the glm module from PyMC3,
which was written by Thomas Wiecki and others, in order to easily specify our
Bayesian linear regression.
Before we Continue let’s describe how MCMC works:
MCMC for Dummies
Let’s start with bayes formula,

P (x|θ)P (θ)
P (θ|x) = (22)
P (x)
P (x) is called the evidence (i.e. the evidence that the data x was generated by
this model). We can compute this quantity by integrating all possible parameter
values:

8
 Z
P (x) = P (x, θ)dθ (23)
Θ
The difficulty here is that for even slightly non-trivial models you can’t just
compute the posterior in a closed-form way.
What about approximate it? For example draw sample from that posterior
we can monte carlo approximate it. Unfortunately to directly sample from that
distribution you not only need to solve the bayes formula, but also invert it,
so that’s even harder (Review how to draw samples from Nomal distribution
(N (µ, σ 2 ))).
Then we might say ”Well, instead let’s construct an ergodic, reversible
Markov chain that has as an equilibrium distribution which matches our poste-
rior distribution”. I’m just kidding, most people wouldn’t say that as it sounds
bat-shit crazy. If you can’t compute it, can’t sample from it, then constructing
that Markov chain with all these properties must be even harder.
The surprising insight though is that this is actually very easy and there
exist a general class of algorithms that do this called Markov chain Monte Carlo
(constructing a Markov chain to do Monte Carlo approximation).
How it works:
Lets generate some data: 20 points from a normal centered around zero.
Our goal will be to estimate the posterior of the mean mu (we’ll assume that
we know the standard deviation to be 1).
Next, we need to define our model. In this case we assume that the data
is normal distribution, i.e. the likelihood of the data is normal. For simplicity,
we assume that we know σ = 1 and we’ll want to infer the posterior of µ. For
each parameter we have to choose a prior. For simplicity, let’s assume a Normal
distribution as a prior for µ. Thus in stats speak our model is:

µ ∼ N ormal(0, 1) (24)

µ|x ∼ N ormal(x; µ, 1) (25)

What is convenient is that for this model, we actually can compute the
posterior analytically. That’s because a normal likelihhod with known std, the
normal prior for µ is conjugate (conjugate here means that our posterior will
follow the same distribution as the prior), so we know that our posterior for µ
is also normal.
However, let’s see how we can sample. At first you find starting parameter
position (can be chosen randomly). Let say that µ current = 1.
Then, you propose to move (jump) from that position somewhere else (that’s
the Markov part). You can be very dumb or very sophisticated about how you
come up with that proposal. The Metropolis sampler is very dumb and just
takes a sample from a normal distribution (no relationship to the normal we
assume for the model) centered around your current mu value (i.e. mucurrent )
with a certain standard deviation (proposalwidth ) that will determine how far
you propose jumps (here we’re use [Link]):

9
 proposal = norm(mucurrent , proposalwidth ).rvs() (26)
Next, you evaluate whether that’s a good place to jump to or not. If the
resulting normal distribution with that muproposed explains the data better than
your old mu, you’ll definitely want to go there. What does ”explains the data
better” mean? We quantify fit by computing the probability of the data, given
the likelihood (normal) with the proposed parameter values (µproposed and a
fixed σ = 1). This can easily be computed by calculating the probability for
each data point using [Link](mu, sigma).pdf (data) and then mul-
tiplying the individual probabilities, i.e. compute the likelihood (usually you
would use log probabilities but we omit this here):
code:

likelihoodcurrent = norm(mucurrent, 1).pdf (data).prod() (27)

likelihoodproposal = norm(muproposal, 1).pdf (data).prod() (28)

compute the priori probability of current and proposed mu.

priorcurrent = norm(mupriormu, mupriorsd).pdf (mucurrent) (29)

priorproposal = norm(mupriormu, mupriorsd).pdf (muproposal) (30)

Up until now, we essentially have a hill-climbing algorithm that would just

propose movements into random directions and only accept a jump if the muproposal
has higher likelihood than mucurrent . Eventually we’ll get to mu = 0 from where
no more moves will be possible. However, we want to get a posterior so we’ll
also have to sometimes accept moves into the other direction. The key trick is
by driving the two probabilities,

paccept = pproposal /pcurrent (31)

We get an acceptance probability. You can already see that if pproposal is larger,
that probability will be > 1 and we will definitely accept. However, if pcurrent
is larger, say twice as large, there’ll be 50% chance of moving there:

accept = [Link]() < paccept (32)

if accept: (33)
curpos = proposal (34)

To summarize the algorithm above:

Given some data {xi }N i=1 and by assuming that data are coming from nor-
mal distribution with known variance σ 2 , our goal is to estimate the posterior
distribution of µ, p(µ|x). Using Bays formula we have:

10
 p(x|µ)p(µ)
p(µ|x) = (35)
p(x)
To estimate this we need to calculate p(x) which can be very challenging in
many use cases. Instead of deriving the posterior. Also if you want to sample
from this distribution, you don’t only need to solve the bayes formula, you need
to also invert it, which also can be hard.
The idea is instead of derive it, we can construct an ergodic, reversible
Markov chain that has an equilibrium distribution that matches our posterior
distribution. There exist a general class of algorithms that do this called Markov
Chain Monte Carlo.
The method is described below:

1. We assume that the data xi are normally distributed:

µ ∼ N (0, 1) (36)
x|µ ∼ N (x; µ, 1) (37)

Actually for this model we can compute the posterior analytically (the
normal prior for µ is conjugate). However, we can still use it an example
to illustrate the MCMC approach.
2. We find starting parameter position (can be chosen randomly). Let’s say
µc = 1
3. We follow the metropolis sampler and we sample from a normal distribu-
tion centered around µc , with a certain std (i.e. σ 1 = 1), mp ∼ N (µc , 1)
4. we check whether the µp explains the data better that the µc . Q To do this we
n
compute the likelihood of the data based on µc , L(x, µc ) = i=1 p(xi |µc )
with
Qn the likelihood of the data based on the newly proposed µp , L(x, µp ) =
i=1 p(xi |µp ). We also compute the prior pdf based on the current µc ,
2 2
f (µc ; µprior , σprior ) and based on the proposed µp , f (µp ; µprior , σprior ).
5. we then calculate the nominator of the Baye’s formula for both µc , pc =
2 2
L(x, µc )·f (µc ; µprior , σprior ) and for µp , pp = L(x, µp )·f (µp ; µprior , σprior )

6. We then calculate:
pp
paccept = (38)
pc
7. We have an acceptance probability. If pp is larger, then paccept > 1 and
we will definitely accept the proposed value µp . However, we might have
pc > pp . So, we accept the proposed µp if paccept > u, where u is a random
sample from uniform distribution over [0, 1).
8. this simple process is repeated with the µc to be the µp .

11
 Taking one step back, note that the above acceptance ratio is the reason this
whole thing works out and we get around the integration. We can show this by
computing the acceptance ratio over normalized posterior and seeing how it’s
equivalent to the acceptance ratio of the unnormalized posterior (let’s say that
µ0 is our current position, and µ is our proposal):
P (x|µ)P (µ)
P (x) P (x|µ)P (µ)
P (x|µ0 )P (µ0 )
= (39)
P (x|µ0 )P (µ0 )
P (x)

In words, dividing the posterior of proposed parameter setting by the pos-

terior of the current parameter setting, P (x) – that nasty quantity we can’t
compute – gets canceled out. So you can intuit that we’re actually dividing the
full posterior at one position by the full posterior at another position (no magic
here). That way, we are visiting regions of high posterior probability relatively
more often than those of low posterior probability.
Now, let’s go back to the Bayesian Linear Regression. We can write the
posterior of the parameters β as

p(y|β, X) · p(β|X)
p(β|y, X) = , (40)
p(y|X)
where p(β|X) = pdf (β|X) is the prior belief of parameter β, p(y|β, X) is the
likelihood of observing y given our prior; and p(y|X) is called the evidence.
We can use a conjugate prior as described above, y ∼ N (Xβ, σ 2 I) and
β ∼ N (β0 , Σβ,0 )

5 Clustering in time series

1. Cluster Analysis on Stock Selection
The scope here is to see whether financial ratios in different dimensions
really add value to stock selection. Also we want to see whether cluster
analysis helps us to get rid of the sea of financial metrics.
First we access the financial ratios and key metrics (e.g. average inventory,
average payables, capex per share, current ratio, debt to assets and etc.)
for the stocks.
We also Download the stock price data. We get the daily price data (close
price). The ratios/metrics in order to capture a company’s performance
or ability from different dimensions are described below:

• Debt to Equity Ratio: Leverage Ratio of a company. It is calculated

by dividing a company’s total liabilities by its shareholder equity
• Dividend Yield: Amount of money a company pays shareholders for
owning a share of its stock dividend by its current stock price.

12
 • Enterprise Value (EV) / EBITDA: A commonly used multiple to
determine the value of a company. It is the earnings before inter-
est taxes depreciation and amortization. A measure of a company;s
overall performance.
• Net Profit Margin - how much net income is generated as a percentage
of revenue
Clustering help us to group similar stocks together based on their char-
acteristics, represented by financial indicators. Here, to perform the clus-
tering we can applied hierarchical clustering, with bottom up approach,
Euclidean Distance and Ward’s method to calculate the similarity.
Hierarchical Clustering Technique: It is one of the popular clustering
techniques. Initially, each data point is considered as an individual cluster.
At each iteration, the similar clusters merge with other clusters until one
cluster or K clusters are formed.

6 Dynamic Time Wrapping

K-means and Euclidean matching algorithm usually fails in time series,
because time series can be of different length, might be related but there
might be a lag between them. The idea to compare arrays with different
length is two build one to many and many to one matches so that the
total distance can be minimized between the two.
There might be series that follow the same pattern, but there might be
of different length. If we apply one-to-one match, shown in the top, the
mapping is not perfectly synced up. DTW overcomes the issue by devel-
oping one-to-many match so that the troughs and peaks with the same
pattern are perfectly matched, and there is no left out for both curves.
In general, DTW is a method that calculated an optimal match between
two given sequences with certain restrictions and rules:
(a) Every index from the first sequence must be matched with the first
index from other sequence and vice versa.
(b) The first index from the first sequence must be matched with the
first index from the other sequence (but it doesn’t have to be its only
match)
(c) The last index from the first sequence must be matched with the last
index from the other sequence (but it does not have to be its only
match)
(d) The mapping of the indeces from the first sequence to the indeces
from the other sequence must be monotonically increasing, and vice
versa. If j > i are indeces from the first sequence, then there must
not be two indices l > k on the other sequence, such that index i is
matched with index l and index j is matched with index k.

13