Please help with fmincon function

1 visualizzazione (ultimi 30 giorni)
dav
dav il 9 Nov 2013
Commentato: dav il 10 Nov 2013
Hi,
I am trying to estimate parameters of an arch model using MLE.
If my code works fine I should get estimates 0.1 and 0.4 ( which I used to create the data)
However, I am having some issues and the estimates I get are wrong.
Can someone help me with this please. Given below are the codes I used.
% data generation
clc;
clear;
p=1;
T = 3000;
ra = zeros(T+2000,1);
seed=123;
rng(seed);
ra = randn(T+2000,1);
epsi=zeros(T+2000,1);
simsig=zeros(T+2000,1);
a0=0.1; a1=0.4;
unvar = a0/(1-a1);
for i = 1:T+2000
if (i==1)
simsig(i) = a0+a1*((a0)/(1-a1));
s=(simsig(i))^0.5;
epsi(i) = ra(i) * s;
else
simsig(i) = a0+ a1*(epsi(i-1))^2;
s=(simsig(i))^0.5;
epsi(i) = ra(i)* s;
end
end
yt = epsi(2001:T+2000);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ESTIMATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
theta0 = [1;1];
A=[0 1];
b=0.99999;
[theta, opt] = fmincon(@(theta)...
lach(theta,yt),theta0,A,b,[],[],0,1);
% function lach
function L = lach(theta,y)
w = theta(1);
alpha = theta(2);
y2 = y.^2;
[T1,K] = size(y2);
ht = zeros(T1,1);
ht(1) = sum(y2)/T1;
for i=2:T1
ht(i)=w + alpha*y2(i-1);
end
sqrtht = sqrt(ht);
x = y./sqrtht;
l = -0.5*log(2*pi) - log(sqrtht) - 0.5*(x.^2);
l=-l;
L = sum(l);
Also, if you can kindly explain the difference between the following two ways of calling fmincon, that will be highly appreciated.
[theta, opt] = fmincon(@(theta)...
lach(theta,yt),theta0,A,b,[],[],0,1);
[theta, opt] = fmincon(@lach,theta0,A,b,[],[],0,1);
thanks
  • When writing these codes I got some help from codes I found online.
  4 Commenti
dav
dav il 10 Nov 2013
thanks
Matt J
Matt J il 10 Nov 2013
Modificato: Matt J il 10 Nov 2013
If my code works fine I should get estimates 0.1 and 0.4 ( which I used to create the data)
What happens when you run fmincon with 0.1 and 0.4 as your initial guess? What is the final value for opt in that case and what was it for
theta =
0.0002
-13.4557

Accedi per commentare.

Risposta accettata

Walter Roberson
Walter Roberson il 10 Nov 2013
[theta, opt] = fmincon(@lach,theta0,A,b,[],[],0,1);
will not have any way of getting y into the routine lach. What you have now,
[theta, opt] = fmincon(@(theta) lach(theta,yt),theta0,A,b,[],[],0,1);
is fine.
The message "Local minimum possible. Constraints satisfied." is a message you would get if it was successful in finding a minimum. fmincon cannot tell whether the minimum it finds is a local minimum or the global minimum so it is being conservative in warning you that what it has find might be a local minima.
The message "fmincon stopped because the size of the current search direction is less than" and so on, is less common. You are getting it because you are hitting against the minimum for the first pararameter. The algorithm has detected that it cannot take two more steps towards the minimum, using the minimum step size that it has been configured to use. You could pass an options structure to specify a smaller minimum step size, but the effect of that would just be to get you ever closer to the boundary 0.
I suspect that you are running into a local minimum. That happens with fmincon(). You can try experimenting with different theta0 (different starting points.) If you had the global optimization toolbox, you could use a routine that would run fmincon() with a variety of different starting points, automatically, looking for the best result out of the set. The global optimization toolbox also has other minimization routines that can be used.
  6 Commenti
Matt J
Matt J il 10 Nov 2013
Modificato: Matt J il 10 Nov 2013
Can you please tell me if i can use fminsearch with constraints ?
There are ways to extend fminsearch to deal with constraints, e.g.,
However, it is not rational that fminsearch does better than fmincon. You might be blaming witches for the crops.
dav
dav il 10 Nov 2013
Thanks a lot for all your help!

Accedi per commentare.

Più risposte (0)

Categorie

Scopri di più su Particle & Nuclear Physics in Help Center e File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by