Skip to content

Instantly share code, notes, and snippets.

@tomaskrehlik
Created August 25, 2014 07:31
Show Gist options
  • Save tomaskrehlik/ce7ac48a523e7b7a0895 to your computer and use it in GitHub Desktop.
Save tomaskrehlik/ce7ac48a523e7b7a0895 to your computer and use it in GitHub Desktop.
Miscellaneous tests for time series.
# Critical values for the Augmented Dickey-Fuller test.
# http://home.cerge-ei.cz/petrz/GDN/crit_values_ADF_KPSS_Perron.pdf
# The p-vals are 0.01, 0.025, 0.05, 0.1
# Number of observations is in the Ts
# Small test only for the df test
# Should print something close to critical value
# using Distributions
# mean([sort([df(cumsum(rand(Normal(), 1000)), "None")[1] for i=1:500])[25] for j=1:200])
using Distributions
none_cv = [ -2.66 -2.26 -1.95 -1.6;
-2.62 -2.25 -1.95 -1.61;
-2.60 -2.24 -1.95 -1.61;
-2.58 -2.23 -1.95 -1.62;
-2.58 -2.23 -1.95 -1.62;
-2.58 -2.23 -1.95 -1.62]
drift_cv = [-3.75 -3.33 -3.00 -2.62;
-3.58 -3.22 -2.93 -2.60;
-3.51 -3.17 -2.89 -2.58;
-3.46 -3.14 -2.88 -2.57;
-3.44 -3.13 -2.87 -2.57;
-3.43 -3.12 -2.86 -2.57]
trend_cv = [-4.38 -3.95 -3.60 -3.24;
-4.15 -3.80 -3.50 -3.18;
-4.04 -3.73 -3.45 -3.15;
-3.99 -3.69 -3.43 -3.13;
-3.98 -3.68 -3.42 -3.13;
-3.96 -3.66 -3.41 -3.12]
Ts = [25, 50, 100, 250, 300]
# The backsolve can be substituted by some reliable OLS machinery when this is available
function adf(ts, typ, lags)
if lags == 0
return adf(ts, typ)
end
Δ = diff(ts)
l = apply(hcat, [Δ[(lags-i+1):(end-i)] for i=1:lags])
if typ=="None"
y = Δ[(lags+1):end]
x = hcat(ts[(lags+1):(end-1)], l)
dfded = 0
cv = none_cv
elseif typ=="Drift"
y = Δ[(lags+1):end]
x = hcat(ts[(lags+1):(end-1)], fill(1.0, length(Δ)-lags-1), l)
dfded = 1
cv = drift_cv
elseif typ=="Trend"
y = Δ[(lags+1):end]
x = hcat(ts[(lags+1):(end-1)], fill(1.0, length(Δ)-lags-1), [1:(length(Δ)-lags-1)], l)
dfded = 2
cv = trend_cv
end
coef = (x'*x) \ x'*y
res = y - x*coef
n = (size(res)[1])
k = lags+1+dfded
σ = sqrt((res'*res)/(n-k))
std = σ[1] * sqrt(inv(x'*x))
DF = (coef ./ diag(std))[1]
AIC = n*log((res'*res)/n) + 2*k
return (DF, AIC)
end
function df(ts, typ)
Δ = diff(ts)
if typ=="None"
y = Δ
x = hcat(ts[1:(end-1)])
dfded = 0
cv = none_cv
elseif typ=="Drift"
y = Δ
x = hcat(ts[1:(end-1)], fill(1.0, length(Δ)-1))
dfded = 1
cv = drift_cv
elseif typ=="Trend"
y = Δ
x = hcat(ts[1:(end-1)], fill(1.0, length(Δ)-1), [1:(length(Δ)-1)])
dfded = 2
cv = trend_cv
end
coef = (x'*x) \ x'*y
res = y - x*coef
n = (size(res)[1])
k = 1+dfded
σ = sqrt((res'*res)/(n-k))
std = σ[1] * sqrt(inv(x'*x))
DF = (coef ./ diag(std))[1]
AIC = n*log((res'*res)/n) + 2*k
return (DF, AIC)
end
function schwertLagMax(ts)
T = size(ts)[1]
return convert(Int64, floor(12*(T/100)^(1/4)))
end
function acov(ts, order)
return cov(ts[(1+order):end], ts[1:(end-order)])
end
# # Small test
# t = rand(Normal(), 100)
# coefs = [0.2, 0.003]
# a = size(coefs)[1]+1
# for i=a:(size(t)[1])
# t[i] = (t[(i-a+1):(i-1)]'*coefs + t[i])[1]
# end
# ljungBox(t, 3)
function ljungBox(ts, order)
n = size(ts)[1]
Q = n*(n+2)*mapreduce(k->acov(ts,k)/(n-k), +, 1:order)
p = ccdf(Chisq(order), Q)
return (Q, p)
end
function boxPierce(ts, order)
n = size(ts)[1]
Q = n*mapreduce(k->acov(ts,k), +, 1:order)
p = ccdf(Chisq(order), Q)
return (Q, p)
end
function durbinWatson(resid)
return sum((resid[2:end]-resid[1:(end-1)])^2)/sum(resid.^2)
end
@milktrader
Copy link

This looks like a great start. Are you interested in doing a PR to HypothesisTests with these? They need a little polishing

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment