1.在matlab ,do maximum likelihood ar1 regression时,是这样子

% use cochrane-orcutt estimates as initial values reso = olsc(y,x);

parm = zeros(nvar+2,1); parm(1,1) = reso.rho; % initial rho parm(2:2+nvar-1,1) = reso.beta; % intiial bhat's parm(2+nvar,1) = reso.sige; % initial sigma

oresult = maxlik('ar1_like',parm,[],y,x);

3.ar1_like(),是用来evaluate ols model with AR1 errors log-likelihood

4.不知是否符合你的需要,需要的话我可提供资料供参考

用matlab最优化函数命令

fminsearch('你的似然函数')；

fminsearch('-（你的似然函数）')；

Below I give an example of GAUSS code estimating garch_t model using "optmum".

Similar for maxlik. Some times I find optmum is more robust than maxlik.¨

but in general maxlik is easier to use.

Note that some part of the code may not be relevant for your purpose.

----------------------

cls; library pgraph,optmum; #include optmum.ext; graphset; data_file2 = "E:\\Timo1\\Stock_index\\SHHChina.xls"; data_range2 = "a302:a1837"; data = SpreadsheetReadM(data_file2,data_range2,1); data1 = lag1(data); data = data[2:rows(data),.]; data1 = data1[2:rows(data1),.]; R = ln(data./data1); y = R-meanc(R);

T_simu_Setting = 1300; // Number of realizations for each simulation Burnin_Setting = 300; // Burn-in Draw_Setting = 500; // Number of Simulations Trim_or_Not = 0; // 1 for trimming, 0 for Not Trimming file_name_suffix = "_500_Untrimmed_fixed_nu.xls"; File_Ex_Plug_in = "E:\\Timo1\\Confi_Regi\\GARCH_t\\GARCH_t_Plug_in_KR_ACF" $+ file_name_suffix; File_Ex_Para_Estimation_Result = "E:\\Timo1\\Confi_Regi\\GARCH_t\\GARCH_t_para_Estimation_Result" $+ file_name_suffix; File_Ex_para_hat_Isoquant = "E:\\Timo1\\Confi_Regi\\GARCH_t\\GARCH_t_para_hat_Isoquant" $+ file_name_suffix; File_Ex_nonparametric_KR_ACF = "E:\\Timo1\\Confi_Regi\\GARCH_t\\GARCH_t_Nonparametric_KR_ACF" $+ file_name_suffix; File_Ex_Esti_Simu = "E:\\Timo1\\Confi_Regi\\GARCH_t\\GARCH_t_Esti_Simu" $+ file_name_suffix; ?;"File_Ex_Plug_in " File_Ex_Plug_in; ?;"File_Ex_para_hat_Isoquant " File_Ex_para_hat_Isoquant; // I found that I need to run this set of codes again and again, so I moved the above command to the // beginning of the program. This will help prevent me from forgetting change some commands when re-run // program with different settings. 2005,11,06

if Trim_or_Not ==1; " The use of the follow part of code is interesting when this part is commented, "; " We use the original log returns possibly with some ourliers(extreme observations). "; " When the simulation cloud can not capture the empirical measures (becasue of the "; " outliers), we then de-comment this code part, and only analyze the dataset with outliers"; " deleted. It is very likely that with the trimmed dataset, the simulation cloud includes "; " the empirical measure of acf and kurtosis very well"; LowQ=0.005; UpperQ=0.995; e_q=LowQ|UpperQ; y_q=quantile(y,e_q); " quantile returns a column vector "; W=y; Q1=(W.>=Y_q[2,.]); Q2=(W.<=Y_q[1,.]); Q=Q1+Q2; " those elements less than the lower or greater than the upper quantile are set as 1"; Q=1-Q; " it becomes zero "; W=W.*Q; " such elements become 0's "; W=miss(W,0); " 0's become missing values "; y=packr(W); " all missing values are deleted "; endif;

LogReturns = y; // this Log Return we be used in the simulation part of this code

b10 = 1.636e-5; b20 = 0.2; b30 = 0.7; b40 = 8; b0 = b10|b20|b30|b40; Kjia= 4; /* we have 4 parameters to be estimated */ _opalgr= 2; /* _opalgr = 1 Steepest decent, default ** 2 BFGS ** 3 Scaled BFGS ** 4 Self_Scaling DFP ** 5 Newton-Raphson ** 6 polak-Ribiere Conjugate Gradient --*/ _opstep= 3; /* _opstep = 1 steplength=1 ** 2 STEPBT, default ** 3 golden steplength ** 4 Brent -------*/ _opshess = 1; Flag72=1; /* 1 for returning objective function for minimization */ {best72,l,s,retcode} = optmum(&GARCH_N,b0); T = rows(y); Flag72 = 2; /* if 2, return the vector individual log likelihood */ /* this is needed to compute the h(theta0,Y_t), the score */ /* in Hamilton 1994 pp.389 */ h0 = gradp(&GARCH_N,best72); /*** reutrn a T*k gradient vector ****/ /*** k is the number of parameter ****/ J_op = zeros(Kjia,Kjia); for i (1,T,1); J_op=J_op+h0[i,.]'*h0[i,.]; endfor; J_op = J_op/T; /* Now J_op is just the J_op in eqn. [13.4.9] in */ /* Hamilton 1994, pp. 389 */ Flag72 = 3; /*-- if 3 return the log-likelihood --*/ hess0 = hessp(&GARCH_N,best72); J_2d = -hess0/T; /* J_2d is the J_2D in Hamilton 1994, pp. 389 [13.4.8] */ cov = inv(J_2d*inv(J_op)*J_2d); cov = cov/T; covG = inv(J_2d)/T; /* This is to create the var-cov matrix using only the */ /* Information Matrix */ Omega = exp(best72[1]); /* revert back to the natural para. for Omega */ alpha = best72[2]^2; /* revert back to the natural para. for Alpha */ w3 = best72[3]; beta = exp(w3)/(1+exp(w3)); /* revert the natural para. Esti. for beat */ nu = exp(best72[4]);

d1 = exp(best72[1]); /* for Delta Method */ d2 = 2*best72[2]; /* for Delta Method */ d3 = exp(w3)/((1+exp(w3))^2); /* for Delta Method */ d4 = exp(best72[4]); /* for Delta Method */

d72 = d1|d2|d3|d4; /* for Delta Method */ d72 = diagrv(zeros(Kjia,Kjia),d72); cov = d72*cov*d72; /* The Delta Method */ covG = d72*covG*d72;?;?;?;?;?;?; est_1= Omega|alpha|beta|nu; cov_1= cov; cov_1_three = cov[1:3,1:3]; // later on we will fix the nu at some rounded off number, so // when we do the simulation, we only draw random number // by viewing the nu_hat as a fixed value, not random. // est_1 and cov-1 are for simulation part ?; " Estimation Results with Robust standard errors."; ?; " Starter estimate Std.Err. t"; "Omega" b10~Omega~sqrt(cov[1,1])~Omega/sqrt(cov[1,1]); "Alpha" b20~alpha~sqrt(cov[2,2])~alpha/sqrt(cov[2,2]); "beta " b30~beta~sqrt(cov[3,3])~beta/sqrt(cov[3,3]); "nu " b40~nu~sqrt(cov[4,4])~nu/sqrt(cov[4,4]);?; "var-cov matrix" cov; ?; "Correlation matrix:" corrvc(cov); ?; " Estimation Results with standard errors from information matrix "; ?; " Starter estimate Std. Error t"; "Omega" b10~Omega~sqrt(covG[1,1])~Omega/sqrt(covG[1,1]); "Alpha" b20~alpha~sqrt(covG[2,2])~alpha/sqrt(covG[2,2]); "beta " b30~beta~sqrt(covG[3,3])~beta/sqrt(covG[3,3]); "nu " b40~nu~sqrt(covG[4,4])~nu/sqrt(covG[4,4]);?; "var-cov matrix " covG;?; "alpha + Beta = " alpha+beta; //-------------------------------------------

proc Garch_N(b); local T,lhd,home,like,i,h,e,alp0,alpha,beta,nu; alp0 = exp(b[1]); /* alp0 > 0 */ alpha = b[2]*b[2]; /* alpha >=0 */ beta = exp(b[3])/(1+exp(b[3])); /* beta is restricted on the interval (0,1)*/ nu = exp(b[4]); /* degree of freedom parameter > 0 */

T = rows(y); h = zeros(T,1); lhd = zeros(T,1); e = y; /* error term in the mean eqn */ h[1,1] = e'*e/T;

/* Initial value of conditional variance is set at */ /* the sample variance of the error term in the */ /* mean equation */ /* h[1,1]=alp0/(1-beta); */ /* as an alternative, the initial value for the */ /* conditional variance can be set at the unconditional */ /* varaicne */ i=2; do while i<=T; h[i,1] = alp0+alpha*e[i-1,1]^2+beta*h[i-1,1]; lhd[i,1] = ln(gamma((nu+1)/2))-ln(gamma(nu/2))-0.5*ln(pi*nu)-0.5*ln(h[i,1])-0.5*(nu+1)*ln(1+e[i,1]^2/(nu*h[i,1])); i=i+1; endo; /* lhd is the vector of sample log-likehood */ like = sumc(lhd); /* sumc becomes the sample log-likelihood, a scalar */ if Flag72 == 1; HOME = -like; /* for minimization with optmum */ elseif Flag72 == 2; HOME = lhd; /* for score function computation */ elseif Flag72 == 3; HOME = like; /* return the sample log-likelihood function, a scalar */ endif; retp(home); endp;

[此贴子已经被作者于2005-12-22 1:44:00编辑过]

感谢各位的热情帮助，可是有没有用stata的人啊？matlab和Gauss我都不会用，只会用sas和stata。另外我想用的似然函数需要模拟，不知道能不能用stata编程完成。再次感谢各位的热情帮助。好温暖啊！（菜鸟发自内心的感激！）