/*to search for maximum in MP problem*/
new;
library optmum, pgraph;

clear d2, X, y;
/* proc calculates rss, y on X with col 26 = logit(d2[.,27])*/
proc rssc(theta);							
	local mmu, msig, z, e;
	
	mmu = theta[1];							@extract mu, sigma@
	msig = theta[2];
	z = 1./(1+exp(-(d2[.,27]-mmu)/msig));	@create logit(d2[.,27]); @ 
	X[.,26] = z;							@replace col 26 of X     @				
	e = y - X*(y/X);						@calculate residuals     @
	retp( 0.5*rows(X)*ln(e'e ) );			@ -ln Likelihood + const @
endp;


load d[]=c:\gausswin\x\BoE\dati1099.csv;	@load Italian data		 @
kT=rows(d);									@number of data points	 @
kT "221*44=" 221*43;						@check load correct	 @
d2 = reshape(d,221,43);						@reshape into data matrix@
"first row";								@print first and last rows@
d2[1,.];
"last row";
d2[221,.];

d2=d2[13:221,.];							@set sample period		@
y = d2[.,1];								@extract dependent y    @
n = rows(d2);								@ no. obs				@
mmu = 0.12; msig=0.21;						@initial values mu, sig @
z = 1./(1+exp(-(d2[.,27]-mmu)/msig));		@logit d2[.,27]			@
"mmu " mmu "msig " msig;"z[1] " z[1];		
X = ones(n,1)~d2[.,2:25]~z~d2[.,30:43];		@initial X@

output file = exoptmum.out reset;			@save O/P in file		@

{theta,f,g,retcode}=optprt(optmum(&rssc, mmu|msig));
											@minimise -ln L, print	@
mmu = theta[1];								@extract mu, sig		@
msig = theta[2];	
z = 1./(1+exp(-(d2[.,27]-mmu)/msig));		@logit d2[.,27]			@
_ptek = "Z.tkf";							@filename to save graph	@
title("Z");									@title for graph		@
xy(seqa(1,1,n),z);							@graph logit			@

varm = inv(_opfhess);						@Var matrix, inv hessian@
se = sqrt(diag(varm));						@ s.e.					@
t = theta./se;								@ t values				@
format 14,6;
"    coeff           se               t";	@ heading				@
theta~se~t;									@coeff~se~t printed		@

output off;