function m = test_opt()
% simulation, estimation, testing

n = 1000;

x1 = normrnd(0,1,n,1);
x2 = unifrnd(0,1,n,1);
e = normrnd(0,1,n,1);
c = 0.5; sigma = 0.2;
b1 = 0.3; b2 = 1-b1;
theta_true = [c b1 b2 sigma]';
y = c + x1*b1 + x2*b2 + sigma*e;

% unrestricted estimation
theta0 = [0.4 0.2 0.8 0.2]';
[theta,FVAL,EXITFLAG,OUTPUT,G,H] = fminunc(@ls,theta0,[],y,[ones(n,1) x1 x2]);
se = sqrt(diag(inv(H)));
disp('Unrestricted estimation')
disp('Estimate   SD')
disp([theta_true theta se])

% restricted estimation
Aeq = [0 1 1 0];
Beq = 1;
ub = [Inf Inf Inf Inf];
lb = [-Inf -Inf -Inf 0];
[theta,logl,exitflag,output,lambda,G,H] = fmincon(@ls,theta0,[],[],Aeq,Beq,lb,ub,[],[],y,[ones(n,1) x1 x2]);
se = sqrt(diag(inv(H)));
disp('Restricted estimation')
disp('Estimate   SD')
disp([theta_true theta se])
m = 1;


function logl = ls(theta,y,X)
[n,p] = size(X);
b = theta(1:p);
sigma = theta(p+1);
r = y-X*b;
logl = -log(2*pi)*n/2 - n*log(sigma) -sum(r.^2)/(2*sigma^2);
logl = -logl;
if isnan(logl) || isinf(logl)
    logl = 1e10;
end



