@ @ @ Routines to estimate the commodity storage model @ @ by maximum likelihood. This file contains all the @ @ routines. Estimation using maxlik @ @ - cspline - Cubic spline interpolation @ @ - grahes - Computes numerical gradient and hessian @ @ - newmax - Maximizes the likelihood using Newton al@ @ gorithm. @ @ - fi - Routine to solve the initial period @ @ - fp - Routine to solve the other periods @ @ - fprice - Initial routine @ @ @ @ @ @ Procedure to interpolate a price @ @ function using cubic splines @ @ @ proc csp(x,y,yp1,ypn); @ @ @ Inputs: @ @ x = values (nodes) @ @ y = f(x) values @ @ yp1 = derivative at the first node @ @ ypn = derivative at the last node @ @ xe = value of x to evaluate @ @ dd = 0 returns the interpolated val @ @ 1 returns the first deriv. @ @ @ local n,y2,u,p,sig,qn,un,a,b,h,i,k,klo,khi,de,ye; n=rows(x); y2=zeros(n,1); u=zeros(n,1); if yp1>1.0e+30; y2[1]=0; u[1]=0; else; y2[1]=-0.5; u[1]=(3./(x[2]-x[1])).*((y[2]-y[1])./(x[2]-x[1])-yp1); endif; i=2; do while i<=n-1; sig=(x[i]-x[i-1])./(x[i+1]-x[i-1]); p=sig.*y2[i-1]+2; y2[i]=(sig-1)./p; u[i]=(6.*((y[i+1]-y[i])./(x[i+1]-x[i])-(y[i]-y[i-1])./(x[i]-x[i-1]))./(x[i+1]- x[i-1])-sig.*u[i-1])./p; i=i+1; endo; if ypn>=1.0e+30; qn=0; un=0; else; qn=0.5; un=(3./(x[n]-x[n-1])).*(ypn-(y[n]-y[n-1])./(x[n]-x[n-1])); endif; y2[n]=(un-qn.*u[n-1])./(qn.*y2[n-1]+1); k=n-1; @ @ @ Generates second derivative @ @ @ do while k>=1; y2[k]=y2[k].*y2[k+1]+u[k]; k=k-1; endo; retp(y2); endp; @ @ @ Procedure to evaluate the x @ @ @ proc csp1(x,y,y2,xe,dd); local klo,khi,k,h,a,b,ye,de,n; n=rows(x); klo=1; khi=n; do while khi-klo>1; k=(khi+klo)./2; if x[k]>xe; khi=k; else; klo=k; endif; endo; h=x[khi]-x[klo]; a=(x[khi]-xe)./h; b=(xe-x[klo])./h; ye=a.*y[klo]+b.*y[khi]+((a^3-a).*y2[klo]+(b^3-b).*y2[khi]).*(h^2)./6; de=(y[khi]-y[klo])./h+((3.*b^2-1).*y2[khi]-(3.*a^2-1).*y2[klo]).*h./6; if dd==0; retp(ye); else; retp(de); endif; endp; @ @ @ Computes the first iteration @ @ @ proc fi(&f1,n,x,d,r,a,b,z); local f1:proc,i,ff,p1,x1,x2,y,i1; if n==0; if a+b*x>=0; y=a+b*x; else; y=0; endif; else; i=1; ff=0; do while i<=rows(z); p1=f1(&f1,n-1,x,d,r,a,b,z); x1=(p1-a)./b; i1=(1-d).*(x-x1); x2=z[i]+i1; ff=ff+(1./rows(z)).*((1-d)./(1-r)).*f1(&f1,n-1,x2,d,r,a,b,z); i=i+1; endo; if ff>a+b.*x; y=ff; else; y=a+b.*x; endif; endif; retp(y); endp; @ @ @ Computes the price function using cubic splines @ @ and functional iteration @ @ @ proc (2) = fp(&f1,&f2,&f3,xa,xn,h,a,b,d,r,z,tolera); @ @ @ Parameters @ @ @ @ f1 = fi.fcp @ @ f2 = csp.fcp @ @ f3 = csp1.fcp @ @ xa = initial node @ @ xn = final node @ @ h = number of nodes @ @ a = demand intercept @ @ b = demand slope @ @ d = spoilage rate @ @ r = real interest rate @ @ z = harvest points (each has same prob.) @ @ tolera = convergence criteria @ local f1:proc,f2:proc,f3:proc,xnodes,fnodes,ftnodes,i,j,fiter,n,ye,xi1,yi, cerro,y2d; @ @ @ Nodes computation @ @ @ xnodes=zeros(h,1); i=1; do while i<=h; xnodes[i,1]=xa+(i-1).*(xn-xa)./(h-1); i=i+1; endo; @ @ @ Initial functional values @ @ @ fnodes=zeros(h,1); n=1; @ number of initial functional iterations @ i=1; do while i<=h; fnodes[i,1]=f1(&f1,n,xnodes[i,1],d,r,a,b,z); i=i+1; endo; y2d=f2(xnodes,fnodes,b,-0.0001); fiter=1; do while fiter==1; ftnodes=zeros(h,1); @ to use in transition nodes @ i=1; do while i<=h; j=1; ye=0; do while j<=rows(z); xi1=z[j]+(1-d).*(xnodes[i,1]-(fnodes[i,1]-a)./b); yi=f3(xnodes,fnodes,y2d,xi1,0); ye=ye+(1./rows(z)).*((1-d)./(1+r)).*yi; j=j+1; endo; if ye>a+b.*xnodes[i,1]; ftnodes[i]=ye; else; ftnodes[i]=a+b.*xnodes[i,1]; endif; i=i+1; endo; @ @ @ Compare ftnodes with fnodes @ @ @ cerro=sumc(abs(fnodes-ftnodes)); "Error=";;format /rd 12,9;cerro; if cerro<=tolera; fiter=0; else; fnodes=ftnodes; y2d=f2(xnodes,fnodes,b,-0.0001); endif; endo; fnodes=ftnodes; retp(xnodes,fnodes); endp; @ @ @ Function definition @ @ @ proc deasol(par,p); @ @ @ Program to solve the competitive storage problem @ @ using the polynomial technique such as in Williams @ @ and Wright (1991). Inelastic supply, and linear @ @ demand. Uses Deaton and Laroque problem @ @ @ local d1,d2,d3,r,tolera,xa,xn,e,pr,h,m,n,s,pe,pp,qq,i,j,kk,streg,finis,finis1, sn0,sn1,bn,rsq,iter,iterc,f0,ll,der,mats,pp1,pp2,qq1,qq2,derq1,derqn, derp1,derpn,tt,ye,va,w,la,ko,ssa,ye1,ye2,ye3,ye4,ye5,d,ny,y2d1,y2d2; @ @ @ Parameters @ @ @ d1=exp(par[1]); d2=-exp(par[2]); r=0.05; d=-r+exp(par[3]); xa=-6; xn=20; h=10; e={-1.755,-1.045,-0.677,-0.386,-0.126,0.126,0.386,0.677,1.045,1.755}; tolera=1.0e-8; {qq,pp}=fp(&fi,&csp,&csp1,xa,xn,h,d1,d2,d,r,e,tolera); m=rows(e); n=rows(pp); @ @ @ Construction of the likelihood @ @ @ mats=sortc(pp~qq,1); pp1=mats[.,1]; qq1=mats[.,2]; mats=sortc(pp~qq,2); pp2=mats[.,1]; qq2=mats[.,2]; derq1=(qq1[2]-qq1[1])./(pp1[2]-pp1[1]); derqn=(qq1[n]-qq1[n-1])./(pp1[n]-pp1[n-1]); derp1=(pp2[2]-pp2[1])./(qq2[2]-qq2[1]); derpn=(pp2[n]-pp2[n-1])./(qq2[n]-qq2[n-1]); y2d1=csp(pp1,qq1,derq1,derqn); y2d2=csp(qq2,pp2,derp1,derpn); tt=rows(p); ll=zeros(tt-1,1); i=1; do while i<=tt-1; ye1=csp1(pp1,qq1,y2d1,p[i],0); ye2=((p[i]-d1)./d2); ye3=0; j=1; @ @ @ Computation of the mean @ @ @ do while j<=m; ny=(1-d).*(ye1-ye2)+e[j]; ye3=ye3+(1./m).*csp1(qq2,pp2,y2d2,ny,0); j=j+1; endo; @ @ @ Computation of the var @ @ @ ye4=0; j=1; do while j<=m; ny=(1-d).*(ye1-ye2)+e[j]; ye5=csp1(qq2,pp2,y2d2,ny,0); ye4=ye4+(1./m).*(ye5^2)-(1./m).*ye3^2; j=j+1; endo; @ @ @ Creation of the log-likelihood @ @ (-) because of Newton @ @ @ ll[i]=-(ln(2.*pi)+ln(ye4)+((p[i+1]-ye3).^2)./(ye4)); i=i+1; endo; retp(ll); endp; @ @ @ Maximimation @ @ @ library maxlik; #include maxlik.ext; maxset; @ @ @ Load prices @ @ @ load pri[88,13]=deadat.txt; @ @ @ Create prices @ @ @ cof=pri[.,1]; whe=pri[.,2]; sug=pri[.,3]; cot=pri[.,4]; ban=pri[.,5]; coc=pri[.,6]; cop=pri[.,7]; jut=pri[.,8]; mai=pri[.,9]; pal=pri[.,10]; ric=pri[.,11]; tea=pri[.,12]; tin=pri[.,13]; p=cof; @ @ @ Initial parameters @ @ @ par0=zeros(3,1); par0[1,1]=-1.332947512; par0[2,1]=-1.847458286; par0[3,1]=-1.66028419; __title="Deaton and Laroque estimation"; _max_Options={bfgs stepbt forward hetcon}; /* _max_Algorithm=5; _max_covpar=1; _max_Linesearch=5; */ {par,l,gradi,hessi,retcode}=maxlik(p,0,&deasol,par0); call maxprt(par,l,gradi,hessi,retcode);