Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Topic: MLE biivariate multiple regression
Replies: 2   Last Post: Oct 30, 2013 4:06 AM

 Messages: [ Previous | Next ]
 Alessandro Beretta Posts: 3 Registered: 10/28/13
Re: MLE biivariate multiple regression
Posted: Oct 30, 2013 3:58 AM

function [coefficients,nLogL,exitflag] = vecar3(Y,Lags,isConstant)

if isConstant
X = [ones(length(Y),1) lagmatrix(Y, Lags)];
else
X = lagmatrix(Y, Lags);
end

Last = Lags(end);
X=X(Last+1:end,:);
Y=Y(Last+1:end,:);

b=inv(X'*X)*X'*Y;

residuals = Y - X * b;

T=length(Y);

cov_res=residuals'*residuals/T;

parameters=[b ; cov_res];

options = optimset('fmincon');
options = optimset(options, 'Display' , 'iter');
options = optimset(options, 'Diagnostics', 'on');
options = optimset(options, 'Algorithm' , 'sqp');
options = optimset(options, 'TolCon' , 1e-7);
tolerance = 2 * optimget(options, 'TolCon');

F= @(parameters) nLogLike(parameters, Y, X, T);

conStruct=linearConstraints(Lags,isConstant,tolerance);

conStruct.nonLinear= @(parameters) nonLinearConstraints(parameters,Lags,tolerance);

[coefficients, nLogL, exitflag]=fmincon(F,parameters, conStruct.A, conStruct.b, conStruct.Aeq,...
conStruct.beq, conStruct.lb, conStruct.ub, conStruct.nonLinear, options);
end

function [nLogL] = nLogLike(parameters, Y, X, T)

beta=parameters(1:end-2,:);
cov=parameters(end-1:end,:);

S=0;
for i=1:T
for k=1:2
for g=1:2
S=S+(Y(i,k)-X(i,:)*beta(:,k))*inv(cov(k,g))*(Y(i,g)-X(i,:)*beta(:,g));
end
end
end

nLogL=0.5*T*2*log(2*pi)+0.5*T*log(det(cov))+0.5*S;
end

function constraints = linearConstraints(Lags,isConstant,tolerance)

nAR=numel(Lags);

if isConstant

LB= [repmat(-10,1,2); repmat(-1+tolerance,nAR*2,2); [tolerance -100;-100 tolerance] ]; % is correct -100 for the LB of the var-cov matrix?
UB= [repmat( 10,1,2); repmat( 1+tolerance,nAR*2,2); repmat( 100,2,2) ]; % and 100 for the UB?

iCount=2:3;
if numel(Lags)>0
degree=Lags(end);
for L=Lags
if (L==1) && (degree>1)
LB(iCount,:)= -2 + tolerance;
UB(iCount,:)= 2 + tolerance;
elseif (L==2) && (degree>2)
LB(iCount,:)= -1.5 + tolerance;
UB(iCount,:)= 1.5 + tolerance;
end
iCount=iCount+2;
end
end

else

LB= [repmat(-1+tolerance,nAR*2,2); [tolerance -100;-100 tolerance] ]; % is correct -100 for the LB of the cov terms?
UB= [repmat( 1+tolerance,nAR*2,2); repmat( 100,2,2) ]; % and 100 for the UB?

iCount=1:2;
if numel(Lags)>0
degree=Lags(end);
for L=Lags
if (L==1) && (degree>1)
LB(iCount,:)= -2 + tolerance;
UB(iCount,:)= 2 + tolerance;
elseif (L==2) && (degree>2)
LB(iCount,:)= -1.5 + tolerance;
UB(iCount,:)= 1.5 + tolerance;
end
iCount=iCount+2;
end
end

end

Aeq=[];
beq=[];

constraints.lb=LB;
constraints.ub=UB;
constraints.A=[];
constraints.b=[];
constraints.Aeq=Aeq;
constraints.beq=beq;
end

function [c,ceq] = nonLinearConstraints(X,Lags,tolerance)

if isConstant
compan=zeros(max(Lags)*2,max(Lags)*2);
indices=2:3;
for i=Lags
compan(1:2,(i*2-1):(i*2))=X(indices,1:2);
indices=indices+2;
end
compan(3:end,1:(end-2))=eye(max(Lags)*2-2);
else
compan=zeros(max(Lags)*2,max(Lags)*2);
indices=1:2;
for i=Lags
compan(1:2,(i*2-1):(i*2))=X(indices,1:2);
indices=indices+2;
end
compan(3:end,1:(end-2))=eye(max(Lags)*2-2);
end

c = abs(eig(compan)) - (1 - tolerance);
ceq=[];
end

Date Subject Author
10/28/13 Alessandro Beretta
10/30/13 Alessandro Beretta
10/30/13 Alessandro Beretta