Andrew Lynch <[email protected]> is using -optimize()- fit a model:
> I have run into a bit of a dilemma with programing a Kalman filter with the
> optimize() procedure. I have the full documentation and have not been able
> to clarify where exactly the optimize() command looks when it optimizes. I
> am minimizing the process below. My goal is to minimize the Log Likelihood
> function, but while I have observed data yt and xt, a and at are unobserved
> and my goal is to generate them. To do so, I need Stata to rerun the do
> loop every time it updates one of the four parameters (b1, b2, q and h). I
> say all that to ask the question, "Does optimize() only look to the v=
> equation to optimize with the parameters or does it look to the entire
> evaluation and move through it in some ordered fashion (thereby honoring the
> do loop)?" I have been unable to get this to run, and I am unsure whether
> it is my dumb coding error at some point or just that optimize() won't do
> what I want it to do.
>
> Any enlightenment would be greatly appreciated!
-optimize()- will call the evaluator function each time it changes the
parameter values. The function value returned by the evaluator is used to
compute numerical derivatives (in the case of type d0 and v0 evaluators) in
addition to the Newton-Raphson step and determining convergence.
Unfortunately, Andrew's evaluator function is referring to undefined objects.
For example, the -n- object (which appears to be a -real scalar-) is not
defined within -maxlik()-. It appears that -n- is probably the number of
elements in -xt- and -yt-; say
n = rows(xt)
The -do- loop in -maxlik()- is building intermediate calculations that
ultimately go into -L-, which appears to be a -real vector-. -L- must be
initialized before it can be "filled-in":
L = J(n,1,.)
Finally, the following vectors need to be initialized before their elements
can be referenced:
at
ft
pt
a
f
p
y
Maybe some of them should be initialized to the column vector of zeros, such
as
at = J(n+1,1,0)
My best suggestion to Andrew is to work with -maxlik()- directly before using
it with -optimize()-. Build a do-file that verifies that -maxlik()- returns
valid -v- values given predetermined parameters values with his -xt-/-yt-
data. Andrew will have to address the issues I mention above, so the
following do-file could be used as a starting point; it runs without a
run-time error, so all Andrew needs to do now is modify -maxlik()- (from
within check.do) until it produces the correct values given -b-, -xt-, and
-yt-.
***** BEGIN: check.do
clear all
// NOTE: I'm just generating some random values here, but you should use
// inputs that will allow you to verify that -mylik()- is producing valid
// results.
set seed 1234
mata:
xt = invnormal(uniform(100,1))
yt = invnormal(uniform(100,1))
// current parameter values
b = J(1,4,0)
void maxlik(todo, b, xt, yt, v, g, h)
{
b1 = b[1]
b2 = b[2]
q = b[3]
// NOTE: this seems like a scale parameter, so you might consider
// treating it as if it were from a log space (as I've done here) to
// ensure that -h- is positive
h = exp(b[4])
n = rows(xt)
L = J(n,1,0)
// NOTE: figure out how the following are supposed to be initialized
at = J(n+1, 1, 0)
ft = J(n+1, 1, 0)
pt = J(n+1, 1, 0)
a = J(n, 1, 0)
f = J(n, 1, 0)
p = J(n, 1, 0)
y = J(n, 1, 0)
i = 1
do {
a[i]=b1*at[i] + b2*xt[i]
p[i]=(b1)^2*pt[i] + q
p[i]=(b1)^2*pt[i] + q
y[i]=a[i]
f[i]=ft[i] + h
at[i+1]=a[i] + ft[i+1]*p[i]*(yt[i]-y[i])
pt[i+1]=p[i] - ft[i+1]*p[i]*p[i]
L[i]=(.5)/f[i] + (.5)*(yt[i]-y[i])*(yt[i]-y[i])/f[i]
i=i+1
} while(i<(n+1))
v=sum(L)
}
// compute the function value
maxlik(., b, xt, yt, v=., ., .)
// display the function value, if you know what the value should be, then
// maybe show the -reldif()- between what you got and what you expect
v
end
***** END: check.do
--Jeff
[email protected]
For reference, here is the Mata code that Andrew submitted:
> /* Optimization Procedure */
>
> void maxlik(todo, b, xt, yt, v, g, h)
> {
> b1=b[1]
> b2=b[2]
> q=b[3]
> h=b[4]
> i=1
>
> do {
> a[i]=b1*at[i] + b2*xt[i]
> p[i]=(b1)^2*pt[i] + q
> y[i]=a[i]
> f[i]=ft[i] + h
> at[i+1]=a[i] + ft[i+1]*p[i]*(yt[i]-y[i])
> pt[i+1]=p[i] - ft[i+1]*p[i]*p[i]
> L[i]=(.5)/f[i] + (.5)*(yt[i]-y[i])*(yt[i]-y[i])/f[i]
> i=i+1
>
> } while(i<(n+1))
>
> v=sum(L)
> }
>
> S = optimize_init()
> optimize_init_evaluator(S, &maxlik())
> optimize_init_evaluatortype(S, "d0")
> optimize_init_params(S, (0,0,0,0))
> optimize_init_which(S, "min")
> optimize_init_argument(S, 1, xt)
> optimize_init_argument(S, 2, yt)
> b=optimize(S)
>
> *
> * For searches and help try:
> * http://www.stata.com/support/faqs/res/findit.html
> * http://www.stata.com/support/statalist/faq
> * http://www.ats.ucla.edu/stat/stata/
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
