mschauer / MicrostructureNoise.jl / 53

Committed 4 Sep 2018 - 11:33 coverage: 98.947% (+0.07%) from 98.876%

Build # 53

Build Type

Pull #4

travis-ci

Committed by

web-flow

Commit Message

Allow observation at start time

Pull Request Pull Request #4: Allow observation at start time

Run Details

33 of 33 new or added lines in 1 file covered. (100.0%)

1 existing line in 1 file now uncovered.

94 of 95 relevant lines covered (98.95%)

426220.91 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

98.95

/src/microstructure.jl


"""
    MicrostructureNoise.Prior(; N, ��1, ��1, ����, ����, ����, ��0, C0)
    MicrostructureNoise.Prior(; kwargs...)

Struct holding prior distribution parameters.
`N` is the number of bins, 
`InverseGamma(��1, ��1)` is the prior of `��[1]` on the first bin,
the prior on the noise variance `��` is `InverseGamma(����, ����)`,
the hidden state ``X_0`` at start time is `Normal(��0, C0)`, 
and `����` is a prior `Distribution` for `��`, 
for example `���� = LogNormal(1., 0.5)`.

Note: All keyword arguments `N, ��1, ��1, ����, ����, ����, ��0, C0`
are mandatory.


Example:
```
prior = MicrostructureNoise.Prior(
N = 40, # number of bins

��1 = 0.0, # prior for the first bin
��1 = 0.0,

���� = 0.3, # noise variance prior InverseGamma(����, ����)
���� = 0.3,

���� = LogNormal(1., 0.5),
��0 = 0.0,
C0 = 5.0
)
```
"""
struct Prior
    N

    ��1
    ��1

    ����
    ����

    ���� 

    ��0
    C0

    function Prior(;��_...)
        �� = Dict(��_)
        return new(
            ��[:N],
            ��[:��1],
            ��[:��1],
            ��[:����],
            ��[:����],
            ��[:����],
            ��[:��0],
            ��[:C0]
            )
    end
end

"""
    MCMC(��::Union{Prior,Dict}, t, y, ��0::Float64, ����, iterations; 
        subinds = 1:1:iterations, ��0::Float64 = 0.0, printiter = 100,
        fixalpha = false, fixeta = false, skipfirst = false) -> td, ��, ��s, ��s, pacc

Run the Markov Chain Monte Carlo procedure for `iterations` iterations,
on data `(t, y)`, where `t` are observation times and `y` are observations.
`��0` is the initial guess for the smoothing parameter `��` (necessary),
`��0` is the initial guess for the noise variance (optional),
and `����` is the stepsize for the random walk proposal for `��`.

Prints verbose output every `printiter` iteration.

Returns `td, ��s, ��s, ��s, pacc`,
`td` is the time grid of the bin boundaries,
`��s`, `��s` are vectors of iterates,
possible subsampled at indices `subinds`,
`��s` is a Matrix with iterates of `��` rows.
`pacc��` is the acceptance probability for the update step of `��`.

`y[i]` is the observation at `t[i]`.

If `skipfirst = true` and `t` and `y` are of equal length,
the observation `y[1]` (corresponding to `t[1]`) is ignored.

If `skipfirst = true` and `length(t) = length(y) + 1`, 
`y[i]` is the observation at `t[i + 1]`.

Keyword args `fixalpha`, `fixeta` when set to `true` allow fixing
`��` and `��` at their initial values. 
"""
function MCMC(��::Union{Prior,Dict}, t, y, ��0::Float64, ����, iterations; subinds = 1:1:iterations, ��0::Float64 = 0.0, printiter = 100, 
    fixalpha = false, fixeta = false, skipfirst = false)
    
    N = ��.N

    ��1 = ��.��1
    ��1 = ��.��1
    ���� = ��.����
    ���� = ��.����
    ���� = ��.���� 
    ��0 = ��.��0
    C0 = ��.C0


    ��x0 = Normal(��0, sqrt(C0))
 
    if skipfirst
        if length(t) == length(y)
            shift = 0
            @info "skip observation y[1] at t[1] (skipfirst == true)"
        elseif length(t) == length(y) + 1
            shift = 1
        else
            throw(DimensionMismatch("expected length(t) or length(t) - 1 observations (skipfirst == true)"))   
        end      
    else
        length(t) != length(y) && throw(DimensionMismatch("expected length(t) observations"))   
        shift = 0
    end

    n = length(t) - 1 # number of increments  
    
    N = ��.N
    m = n �� N
    


    ���� = InverseGamma(��.����, ��.����)


    # Initialization
    if shift == 1
        x = [y[1]; y]
    else
        x = copy(y)
    end
    �� = ��0
    �� = ��0

    �� = zeros(N)
    �� = zeros(N - 1)



    Z = zeros(N)
    ii = Vector(undef, N) # vector of start indices of increments
    td = zeros(N+1)
    for k in 1:N
        if k == N
            ii[k] = 1+(k-1)*m:n # sic!
        else
            ii[k] = 1+(k-1)*m:(k)*m
        end

        tk = t[ii[k]]
        td[k] = tk[1]
        Z[k] = sum((x[i+1] - x[i]).^2 ./ (t[i+1]-t[i]) for i in ii[k])
        ��[k] = mean(InverseGamma(��1 + length(ii[k])/2, ��1 + Z[k]/2))
    end
    td[end] = t[end]

    acc = 0
    ��s = Float64[]
    si = 1
 
    �� = zeros(n+1)
    C = zeros(n+1)
    ��s = Any[]
    ��s = Float64[]

    samples = zeros(N, length(subinds))

    for iter in 1:iterations
        # update Zk (necessary because x changes)
        if !(�� == 0.0 && fixeta)
            for k in 1:N
                Z[k] = sum((x[i+1] - x[i]).^2 ./ (t[i+1]-t[i]) for i in ii[k])
            end
        end

        # sample chain
        for k in 1:N-1
            ��[k] = rand(InverseGamma(�� + ��, (��/��[k] + ��/��[k+1])))
        end
        for k in 2:N-1
            ��[k] = rand(InverseGamma(�� + �� + m/2, (��/��[k-1] + ��/��[k] + Z[k]/2)))
        end
        ��[1] = rand(InverseGamma(��1 + �� + m/2, ��1 + ��/��[1] + Z[1]/2))
        ��[N] = rand(InverseGamma(�� + length(ii[N])/2, ��/��[N-1] + Z[N]/2))
   
        if !fixalpha
            # update parameter alpha using Wilkinson II
            ���� = �� + ����*randn()
            while ���� < eps()
                ���� = �� + ����*randn()
            end

            lq = logpdf(����, ��)
            lq += (2*(N-1))*(��*log(��) - lgamma(��))
            s = sum(log(��[k-1]*��[k]*��[k-1]*��[k-1]) + (1/��[k-1] + 1/��[k])/��[k-1]  for k = 2:N)
            lq += -��*s

            lq�� = logpdf(����, ����)
            lq�� += (2*(N-1))*(����*log(����) - lgamma(����))
            lq�� += -����*s

            mod(iter, printiter) == 0 && print("$iter \t �� ", ����)
            if rand() < exp(lq�� - lq)*cdf(Normal(0, ����), ��)/cdf(Normal(0, ����), ����) # correct for support
                acc = acc + 1
                �� = ����
                mod(iter, printiter) == 0 && print("���")
            end
            

        end
        if !fixeta
            # update eta
            z = sum((x[i] - y[i - shift])^2 for i in (1 + skipfirst):(n+1))
            �� = rand(InverseGamma(���� + n/2, ���� + z/2))

            mod(iter, printiter) == 0 && print("\t �����", ���(��))


        end
        mod(iter, printiter) == 0 && println()

        # Sample x from Kalman smoothing distribution

        # Forward pass
        if �� == 0.0 && fixeta
            # do nothing
        else 
            if skipfirst # ignore observation y[1]
                C[1] = C0
                ��[1] = ��0
                ��i = ��0
                Ci = C0
            else
                wi = 0.0 
                Ki = C0/(C0 + ��)
                ��i =  ��0 + Ki*(y[1] - ��0) # shift = 0
   
                Ci = Ki*��
                C[1] = Ci
                ��[1] = ��i
            end
            for k in 1:N
                iik = ii[k]
                for i in iik # from 1 to n
                    wi = ��[k]*(t[i+1] - t[i])
                    Ki = (Ci + wi)/(Ci + wi + ��) # Ci is still C(i-1)
                    ��i =  ��i + Ki*(y[i + 1 - shift] - ��i)

                    Ci = Ki*��
                    C[i+1] = Ci # C0 is at C[1], Cn is at C[n+1] etc.
                    ��[i+1] = ��i
                end
            end
            hi = ��i
            Hi = Ci

            # Backward pass
            x[end] = rand(Normal(hi, sqrt(Hi)))

            for k in N:-1:1
                iik = ii[k]
                for i in iik[end]:-1:iik[1] # n to 1
                    wi1 = ��[k]*(t[i+1] - t[i])
                    Ci = C[i]
                    ��i = ��[i]


                    Hi = (Ci*wi1)/(Ci + wi1)
                    hi = ��i + Ci/(Ci + wi1)*(x[i+1] - ��i) # x[i+1] was sampled in previous step
                    x[i] = rand(Normal(hi, sqrt(Hi)))
                end
            end
        end
        if iter in subinds
            push!(��s, ��)
            push!(��s, ��)
            samples[:, si] = ��
            si += 1
        end
    end

    td, samples, ��s, ��s, round(acc/iterations, digits=3)
end

"""
```
struct Posterior
    post_t # Time grid of the bins
    post_qlow # Lower boundary of marginal credible band
    post_median # Posterior median
    post_qup # Upper boundary of marginal credible band
    post_mean # Posterior mean of `s^2`
    post_mean_root # Posterior mean of `s`
    qu # `qu*100`-% marginal credible band
end
```

Struct holding posterior information for squared volatility `s^2`.
"""
struct Posterior
    post_t
    post_qlow
    post_median
    post_qup
    post_mean
    post_mean_root
    qu
end

"""
    posterior_volatility(td, ��s; burnin = size(��s, 2)��3, qu = 0.90)

Computes the `qu*100`-% marginal credible band for squared volatility `s^2` from `��`.

Returns `Posterior` object with boundaries of the marginal credible band,
posterior median and mean of `s^2`, as well as posterior mean of `s`.
"""
function posterior_volatility(td, samples; burnin = size(samples, 2)��3, qu = 0.90)
    p = 1.0 - qu 
    A = view(samples, :, burnin:size(samples, 2))
    post_qup = mapslices(v-> quantile(v, 1 - p/2), A, dims=2)
    post_mean = mean(A, dims=2)
    post_mean_root = mean(sqrt.(A), dims=2)
    post_median = median(A, dims=2)
    post_qlow = mapslices(v-> quantile(v,  p/2), A, dims=2)
    Posterior(
        td,
        post_qlow,
        post_median,
        post_qup,
        post_mean,
        post_mean_root,
        qu
    )
end

"""
    piecewise(t, y, [endtime]) -> t, xx

If `(t, y)` is a jump process with piecewise constant paths and jumps 
of size `y[i]-y[i-1]` at `t[i]`, piecewise returns coordinates path 
for plotting purposes. The second argument
allows to choose the right endtime of the last interval.
"""
function piecewise(t_, y, tend = t_[end])
    t = [t_[1]]
    n = length(y)
    append!(t, repeat(t_[2:n], inner=2))
    push!(t, tend)
    t, repeat(y, inner=2)
end






1
2	"""
3	MicrostructureNoise.Prior(; N, ��1, ��1, ��, ��, ��, ��0, C0)
4	MicrostructureNoise.Prior(; kwargs...)
5
6	Struct holding prior distribution parameters.
7	`N` is the number of bins,
8	`InverseGamma(��1, ��1)` is the prior of `��[1]` on the first bin,
9	the prior on the noise variance `��` is `InverseGamma(��, ��)`,
10	the hidden state ``X_0`` at start time is `Normal(��0, C0)`,
11	and `��` is a prior `Distribution` for `��`,
12	for example `�� = LogNormal(1., 0.5)`.
13
14	Note: All keyword arguments `N, ��1, ��1, ��, ��, ��, ��0, C0`
15	are mandatory.
16
17
18	Example:
19	```
20	prior = MicrostructureNoise.Prior(
21	N = 40, # number of bins
22
23	��1 = 0.0, # prior for the first bin
24	��1 = 0.0,
25
26	�� = 0.3, # noise variance prior InverseGamma(��, ��)
27	�� = 0.3,
28
29	�� = LogNormal(1., 0.5),
30	��0 = 0.0,
31	C0 = 5.0
32	)
33	```
34	"""
35	struct Prior
36	N
37
38	��1
39	��1
40
41	��
42	��
43
44	��
45
46	��0
47	C0
48
49	function Prior(;��_...)
50	�� = Dict(��_)	4×
51	return new(	2×
52	��[:N],
53	��[:��1],
54	��[:��1],
55	��[:��],
56	��[:��],
57	��[:��],
58	��[:��0],
59	��[:C0]
60	)
61	end
62	end
63
64	"""
65	MCMC(��::Union{Prior,Dict}, t, y, ��0::Float64, ��, iterations;
66	subinds = 1:1:iterations, ��0::Float64 = 0.0, printiter = 100,
67	fixalpha = false, fixeta = false, skipfirst = false) -> td, ��, ��s, ��s, pacc
68
69	Run the Markov Chain Monte Carlo procedure for `iterations` iterations,
70	on data `(t, y)`, where `t` are observation times and `y` are observations.
71	`��0` is the initial guess for the smoothing parameter `��` (necessary),
72	`��0` is the initial guess for the noise variance (optional),
73	and `��` is the stepsize for the random walk proposal for `��`.
74
75	Prints verbose output every `printiter` iteration.
76
77	Returns `td, ��s, ��s, ��s, pacc`,
78	`td` is the time grid of the bin boundaries,
79	`��s`, `��s` are vectors of iterates,
80	possible subsampled at indices `subinds`,
81	`��s` is a Matrix with iterates of `��` rows.
82	`pacc��` is the acceptance probability for the update step of `��`.
83
84	`y[i]` is the observation at `t[i]`.
85
86	If `skipfirst = true` and `t` and `y` are of equal length,
87	the observation `y[1]` (corresponding to `t[1]`) is ignored.
88
89	If `skipfirst = true` and `length(t) = length(y) + 1`,
90	`y[i]` is the observation at `t[i + 1]`.
91
92	Keyword args `fixalpha`, `fixeta` when set to `true` allow fixing
93	`��` and `��` at their initial values.
94	"""
95	function MCMC(��::Union{Prior,Dict}, t, y, ��0::Float64, ��, iterations; subinds = 1:1:iterations, ��0::Float64 = 0.0, printiter = 100,
96	fixalpha = false, fixeta = false, skipfirst = false)
97
98	N = ��.N	20×
99
100	��1 = ��.��1
101	��1 = ��.��1
102	�� = ��.��
103	�� = ��.��
104	�� = ��.��
105	��0 = ��.��0
106	C0 = ��.C0
107
108
109	��x0 = Normal(��0, sqrt(C0))	12×
110
111	if skipfirst	12×
112	if length(t) == length(y)	8×
113	shift = 0	4×
114	@info "skip observation y[1] at t[1] (skipfirst == true)"	4×
115	elseif length(t) == length(y) + 1	8×
116	shift = 1	8×
117	else
118	throw(DimensionMismatch("expected length(t) or length(t) - 1 observations (skipfirst == true)"))	8×
119	end
120	else
121	length(t) != length(y) && throw(DimensionMismatch("expected length(t) observations"))	16×
122	shift = 0	160×
123	end
124
125	n = length(t) - 1 # number of increments	16×
126
127	N = ��.N	152×
128	m = n �� N	160×
129
130
131
132	�� = InverseGamma(��.��, ��.��)	160×
133
134
135	# Initialization
136	if shift == 1	160×
137	x = [y[1]; y]	160×
138	else
139	x = copy(y)	8×
140	end
141	�� = ��0
142	�� = ��0
143
144	�� = zeros(N)
145	�� = zeros(N - 1)
146
147
148
149	Z = zeros(N)
150	ii = Vector(undef, N) # vector of start indices of increments
151	td = zeros(N+1)
152	for k in 1:N	8×
153	if k == N	8×
154	ii[k] = 1+(k-1)*m:n # sic!	6,120×
155	else
156	ii[k] = 1+(k-1)m:(k)m	6,120×
157	end
158
159	tk = t[ii[k]]	61,200×
160	td[k] = tk[1]	6,120×
161	Z[k] = sum((x[i+1] - x[i]).^2 ./ (t[i+1]-t[i]) for i in ii[k])	116,280×
162	��[k] = mean(InverseGamma(��1 + length(ii[k])/2, ��1 + Z[k]/2))	6,120×
163	end
164	td[end] = t[end]	110,160×
165
166	acc = 0	6,120×
167	��s = Float64[]	3,060×
168	si = 1	3,060×
169
170	�� = zeros(n+1)	3,060×
171	C = zeros(n+1)	6,120×
UNCOV 172	��s = Any[]	!
173	��s = Float64[]	3,060×
174
175	samples = zeros(N, length(subinds))	6,120×
176
177	for iter in 1:iterations	6,120×
178	# update Zk (necessary because x changes)
179	if !(�� == 0.0 && fixeta)	3,060×
180	for k in 1:N	3,060×
181	Z[k] = sum((x[i+1] - x[i]).^2 ./ (t[i+1]-t[i]) for i in ii[k])	6,120×
182	end
183	end
184
185	# sample chain
186	for k in 1:N-1	3,060×
187	��[k] = rand(InverseGamma(�� + ��, (��/��[k] + ��/��[k+1])))	6,150×
188	end
189	for k in 2:N-1	9,180×
190	��[k] = rand(InverseGamma(�� + �� + m/2, (��/��[k-1] + ��/��[k] + Z[k]/2)))
191	end
192	��[1] = rand(InverseGamma(��1 + �� + m/2, ��1 + ��/��[1] + Z[1]/2))
193	��[N] = rand(InverseGamma(�� + length(ii[N])/2, ��/��[N-1] + Z[N]/2))	5,924×
194
195	if !fixalpha	6,120×
196	# update parameter alpha using Wilkinson II
197	�� = �� + ��*randn()	3,060×
198	while �� < eps()	3,060×
199	�� = �� + ��*randn()	6,150×
200	end
201
202	lq = logpdf(��, ��)	9,180×
203	lq += (2(N-1))(��*log(��) - lgamma(��))	3,060×
204	s = sum(log(��[k-1]��[k]��[k-1]*��[k-1]) + (1/��[k-1] + 1/��[k])/��[k-1] for k = 2:N)	3,060×
205	lq += -��*s	60×
206
207	lq�� = logpdf(��, ��)	60×
208	lq�� += (2(N-1))(��*log(��) - lgamma(��))
209	lq�� += -��*s	60×
210
211	mod(iter, printiter) == 0 && print("$iter \t �� ", ��)
212	if rand() < exp(lq�� - lq)*cdf(Normal(0, ��), ��)/cdf(Normal(0, ��), ��) # correct for support	3,000×
213	acc = acc + 1	3,000×
214	�� = ��	3,000×
215	mod(iter, printiter) == 0 && print("��")	3,000×
216	end
217
218
219	end
220	if !fixeta	3,000×
221	# update eta
222	z = sum((x[i] - y[i - shift])^2 for i in (1 + skipfirst):(n+1))	6,120×
223	�� = rand(InverseGamma(�� + n/2, �� + z/2))	61,200×
224
225	mod(iter, printiter) == 0 && print("\t ��", ��(��))	61,200×
226
227
228	end
229	mod(iter, printiter) == 0 && println()	3,060,000×
230
231	# Sample x from Kalman smoothing distribution
232
233	# Forward pass
234	if �� == 0.0 && fixeta	6,120,000×
235	# do nothing
236	else
237	if skipfirst # ignore observation y[1]	3,060,000×
238	C[1] = C0	3,060,000×
239	��[1] = ��0	3,060,000×
240	��i = ��0	3,060,000×
241	Ci = C0	3,060×
242	else
243	wi = 0.0
244	Ki = C0/(C0 + ��)	3,060×
245	��i = ��0 + Ki*(y[1] - ��0) # shift = 0	6,120×
246
247	Ci = Ki*��	61,200×
248	C[1] = Ci	61,200×
249	��[1] = ��i	3,060,000×
250	end
251	for k in 1:N	3,060,000×
252	iik = ii[k]	3,060,000×
253	for i in iik # from 1 to n	3,060,000×
254	wi = ��[k]*(t[i+1] - t[i])	3,060,000×
255	Ki = (Ci + wi)/(Ci + wi + ��) # Ci is still C(i-1)	3,060,000×
256	��i = ��i + Ki*(y[i + 1 - shift] - ��i)	3,060×
257
258	Ci = Ki*��
259	C[i+1] = Ci # C0 is at C[1], Cn is at C[n+1] etc.
260	��[i+1] = ��i	3,060×
261	end
262	end
263	hi = ��i	6,120×
264	Hi = Ci	24×
265
266	# Backward pass
267	x[end] = rand(Normal(hi, sqrt(Hi)))
268
269	for k in N:-1:1
270	iik = ii[k]
271	for i in iik[end]:-1:iik[1] # n to 1
272	wi1 = ��[k]*(t[i+1] - t[i])
273	Ci = C[i]
274	��i = ��[i]
275
276
277	Hi = (Ci*wi1)/(Ci + wi1)
278	hi = ��i + Ci/(Ci + wi1)*(x[i+1] - ��i) # x[i+1] was sampled in previous step
279	x[i] = rand(Normal(hi, sqrt(Hi)))
280	end
281	end
282	end
283	if iter in subinds
284	push!(��s, ��)
285	push!(��s, ��)
286	samples[:, si] = ��
287	si += 1
288	end
289	end
290
291	td, samples, ��s, ��s, round(acc/iterations, digits=3)
292	end
293
294	"""
295	```
296	struct Posterior
297	post_t # Time grid of the bins
298	post_qlow # Lower boundary of marginal credible band
299	post_median # Posterior median
300	post_qup # Upper boundary of marginal credible band
301	post_mean # Posterior mean of `s^2`
302	post_mean_root # Posterior mean of `s`
303	qu # `qu*100`-% marginal credible band
304	end
305	```
306
307	Struct holding posterior information for squared volatility `s^2`.
308	"""
309	struct Posterior
310	post_t	2×
311	post_qlow
312	post_median
313	post_qup
314	post_mean
315	post_mean_root
316	qu
317	end
318
319	"""
320	posterior_volatility(td, ��s; burnin = size(��s, 2)��3, qu = 0.90)
321
322	Computes the `qu*100`-% marginal credible band for squared volatility `s^2` from `��`.
323
324	Returns `Posterior` object with boundaries of the marginal credible band,
325	posterior median and mean of `s^2`, as well as posterior mean of `s`.
326	"""
327	function posterior_volatility(td, samples; burnin = size(samples, 2)��3, qu = 0.90)
328	p = 1.0 - qu	6×
329	A = view(samples, :, burnin:size(samples, 2))
330	post_qup = mapslices(v-> quantile(v, 1 - p/2), A, dims=2)	86×
331	post_mean = mean(A, dims=2)	2×
332	post_mean_root = mean(sqrt.(A), dims=2)	2×
333	post_median = median(A, dims=2)	2×
334	post_qlow = mapslices(v-> quantile(v, p/2), A, dims=2)	86×
335	Posterior(	2×
336	td,
337	post_qlow,
338	post_median,
339	post_qup,
340	post_mean,
341	post_mean_root,
342	qu
343	)
344	end
345
346	"""
347	piecewise(t, y, [endtime]) -> t, xx
348
349	If `(t, y)` is a jump process with piecewise constant paths and jumps
350	of size `y[i]-y[i-1]` at `t[i]`, piecewise returns coordinates path
351	for plotting purposes. The second argument
352	allows to choose the right endtime of the last interval.
353	"""
354	function piecewise(t_, y, tend = t_[end])
355	t = [t_[1]]	6×
356	n = length(y)
357	append!(t, repeat(t_[2:n], inner=2))	2×
358	push!(t, tend)
359	t, repeat(y, inner=2)	4×
360	end
361
362
363
364
365

mschauer / MicrostructureNoise.jl / 53

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous