Miscellaneous tests for time series.

gistfile1.jl

# 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

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