/*******************************/ /* Examples from the help file */ /*******************************/ /* Help example one */ mata: mata clear /* Regular function */ real matrix func(real matrix x,real scalar i) { if (i==1) return(1:-2*x[,1]:*x[,2]:-3:*x[,1]:^2) if (i==2) return(1:-x[,1]:^2:-3:*x[,2]:^2) } /* Regular jacobian */ real matrix jac(real matrix x,real scalar i, real scalar j) { if (i==1 & j==1) return(-2*x[,2]-6*x[,1]) if (i==1 & j==2) return(-2*x[,1]) if (i==2 & j==1) return(-2*x[,1]) if (i==2 & j==2) return(-6*x[,2]) } /* Interval version of function */ real matrix func_I(real matrix x,real scalar i, real scalar d) { real matrix X1, X2, A, I, B, C, D X1=x[,1::2] X2=x[,3::4] if (i==1) { A=2*int_mult(X1,X2,d) I=J(rows(A),2,1) B=int_sub(I,A,d) C=3*int_mult(X1,X1) D=int_sub(B,C,d) return(D) } if (i==2) { A=int_mult(X1,X1) B=3*int_mult(X2,X2) I=J(rows(A),2,1) C=int_sub(I,A,d) D=int_sub(C,B,d) return(D) } } /* Interval version of jacobian */ real matrix jac_I(real matrix x,real scalar i,real scalar j, real scalar d) { real matrix X1, X2, A, B, C real scalar r X1=x[,1::2] X2=x[,3::4] r=rows(X1) if (i==1 & j==1) { A=int_mult(J(r,2,-2),X2,d) B=int_mult(J(r,2,-6),X1,d) C=int_add(A,B,d) return(C) } if (i==1 & j==2) { A=int_mult(J(r,2,-2),X1,d) return(A) } if (i==2 & j==1) { A=int_mult(J(r,2,-2),X1,d) return(A) } if (i==2 & j==2) { A=J(r,2,-6) B=int_mult(A,X2,d) return(B) } } /* Initialize the problem */ Prob=int_prob_init() int_prob_args(Prob,2) int_prob_f_Iform(Prob,&func_I()) int_prob_jac_Iform(Prob,&jac_I()) int_prob_f(Prob,&func()) int_prob_jac(Prob,&jac()) int_prob_ival(Prob,(-100,100) \ (-100,100)) /* Solve the problem */ int_solve(Prob) int_prob_ints_vals(Prob) int_newton_iter(Prob) int_prob_pts_vals(Prob) end /* Extending the example a little bit */ /* Using routines with random start points */ /* And then just using parallel newton iteration */ mata Prob2=int_prob_init() int_prob_args(Prob2,2) int_prob_f(Prob2,&func()) int_prob_jac(Prob2,&jac()) int_prob_ival(Prob2,(-100,100) \ (-100,100)) /* Different steps from this point on */ randpoints=-1:+2*runiform(10,2) /* Random points */ int_prob_init_pts(Prob2,randpoints) int_prob_method(Prob2,"newton") /* Solution */ int_newton_iter(Prob2) int_prob_pts_vals(Prob2) end /***********************************************/ /* Second problem from the help file */ /* Numerical solution of Diamond's (1982, JPE) */ /* Search model of unemployment */ /***********************************************/ mata /* For the required density function, use a pareto distribution. */ /* Hence: */ /* CDF is F(x)=1-x^alpha, a<0, x>1 */ /* PDF is f(x)= -alphax^(alpha-1), a<0, x>1 */ /* Cond. exp is E[x|x