/**************************************************************************** n_s.p Estimación de las curvas Spot y Forward utilizando el método de Nelson y Siegel. Programa original de Lars Svensson Modificado por Luis Fernando Melo Velandia y Diego Vasquez ****************************************************************************/ library pgraph,maxlik; /* Switches */ obs = 1; /* obs in data file to estimate. 0 = run through all */ nsx = 0; /* 1= extended NS */ _yield = 1; /* 0= min price error, 1= min yield error */ restr = 1; /* 1 = restriction, _s0 /= 0 */ y0 = 0; /* 0=> _s0 = o/n rate, else _s0 = 100*ln(1+y0/100) */ quest = 1; /* 1= ask about saving data */ opt = 1; /* 1 = optimize */ plgr1 = 1; /* 1 = plot Yields to maturity etc */ plgr2 = 0; /* 1 = plot Observed and estimated bond prices */ plgr3 = 0; /* 1 = plot Spot and forward rates */ plgr4 = 0; /* 1 = plot Discount function */ prf = 0; /* 1 = print output file */ prgr = 0; /* 1 = print graphs */ xtics1 = 1; /* 0,1(2000-2009) or 2 (2000-2004) */ ytics1 = 1; /* 0,1 or 2 */ diagnost_ml = 0; /* 1 = ML diagnostic */ toler = 1e-4; /* Default 1e-5 , for convergence in ML*/ max_algorit = 0 ; @ 0->Default, 2->BFGLS, 3->DFP, 4->NEWTON; @ max_mlcov = 0; @ 0->from Hessian(Def.) , 2->from QML; @ Discount_ci = 0; @ 1-> Plot discount function with C.I; 0-> Plot discoint function w/o C.I.; @ /******************************************************************************/ grid = 0; /* 1 = grid */ asy = 1; /* 1 = plot asymptot in graph1 */ mmin = 0.001; /* Mimimum maturity plotted */ mmax= 30; /* Maximum maturity plotted */ nm1 = 501; /* Number of maturities plotted */ errbar = 0; /* 1 = errorbars */ errbarp = 0; /* 1 = errorbars in graph 3*/ ebmin = 1/360; /* Minimum maturity of errorbars */ ebstep = 1; /* Minimum maturity of errorbars */ neb = 8; /* Number of error bars */ ebfact1 = .1; ebfact2 = 2; /* Errorbar factors */ q = 0.95; /* Confidence level for confidence interval */ tol = 0.00001; /* Tolerance level for confidence interval */ mess = 0; /* 1 = LS/MB message and date */ pdate = 1; /* 1 = date and time on graph */ _ncomp = 1; /* compounding and coupon frequency, /yr */ _ytol = 0.000001; /* tolerance for YTM2 */ multb = 0; /* 1 = multiply tb observations */ /* DO NOT Use */ simplon = 0; /* 1= o/n simple rate, calculate compounded rate */ simpltb = 0; /* 1= tb simple rates, calculate compounded rates */ output_file =0; /* 1= Generates 3 files: 1. maturity, Estimated spot, forward and discunt funcions. 2. observed and estimated yield to maturity. 3. observed and estimated bond prices. */ /******************************************************************************/ @==============================================================================@ /*************************************************************************** YAC2 Purpose: Compute yields to maturity compounded _ncomp per year from simple rates, if coupon is zero. Format: ytm = YAC2(c,m,y); Input: c Nx1 vector coupons, %/yr m Nx1 vector times to maturity, yrs y Nx1 vector yields to maturity, %/yr simple yields if c = 0; compounded _ncomp per year if c > 0 Output: ytm Nx1 vector annualized yields to maturity, compounded _ncomp per year, %/yr Global: _ncomp scalar coupon and compounding frequency Remark: Inputs are column vectors, not row vectors as in YAC ----------------------------------------------------------------------------*/ proc YAC2(c,m,y); local h,nm,zz,ytm; h = _ncomp; nm = rows(c); zz = zeros(nm,1); ytm = (c .== zz) .* (100*h*( (1+(y/100).*m)^(1./(m*h)) - 1 )) + (c .> zz) .* y; retp( ytm ); endp; @==============================================================================@ /******************************************************************************/ /* Starting values */ /* b1 b2 b3 b4 b5 b6 */ if restr; bst = 15.8| 1| 2| -1.56| 0.58; @bst = 13.586041 | 6.824524| 1.251817| 0| 0.1;@ @bst = 17| 1| 2| 0| 0.1;@ else; bst = 15.8| -4.8| 1| 2| 0 | 0.1 ; endif; /* @Valores Iniciales de Svensson@ if restr; bst = 7.8| 3| 3.4| -1.56| 0.58; else; bst = 7.2| -0.6| -6| 2.3| 0 | 0.1 ; endif; */ /******************************************************************************/ if opt; _cf = 0; _tp = 0; _ytest = 0; _ptest = 0; /* declare globals */ @ #include NS2.SRC ; #include NS2X.SRC; @ endif; /******************************************************************************/ /* Input */ loadm cm0[] = c:\archivos\term_str\svensson\col\COLsemc.d2; cm0 = reshape(cm0,rows(cm0)/14,14); /* COL?.d2 have 14 cols */ loadm mm0[] = c:\archivos\term_str\svensson\col\COLsemm.d2; mm0 = reshape(mm0,rows(mm0)/14,14); /* COL?.d2 have 14 cols */ loadm ym0[] = c:\archivos\term_str\svensson\col\COLsemy.d2; ym0 = reshape(ym0,rows(ym0)/14,14); /* COL?.d2 have 14 cols */ nobs = rows(ym0); cm = cm0[obs,.]; mm = mm0[obs,.]; ym = ym0[obs,.]; /* Ejercicio sin O/N */ @cm = cm0[obs,1 3:14];@ @mm = mm0[obs,1 3:14];@ @ym = ym0[obs,1 3:14];@ /******************************************************************************/ tstart = date; output file = c:\archivos\term_str\svensson\col\COL.out reset; /******************************************************************************/ /* Adjust for earlier convention that o/n rate has maturity 0 */ if mm[2] == 0; mm[2] = mm[1] + 1/365; endif; tdate = mm[1]; /* Trade date */ ttm = gettime(tdate); /* Trade time */ md = mm[2:cols(mm)]; md = packr(md')'; nm = cols(md); tm = gettime(md); tm = tm - ttm; c = cm[2:cols(cm)]; c = packr(c')'; y = ym[2:cols(ym)]; y = packr(y')'; /* Change simple rates to annualized rates, _ncomp frequency */ if simplon; y[1,1] = YAC2(c[1,1],tm[1,1],y[1,1]); endif; if simpltb; y[1,2:nm] = YAC2(c[1,2:nm]',tm[1,2:nm]',y[1,2:nm]')'; endif; /*** Mulb borrado **/ _s0 = 0; /* Global for NEST2 proc */ if restr; _s0 = 100*_ncomp*ln(1+y[1]/(100*_ncomp)); /* global for NEST2 proc, continuously comp */ if y0; _s0 = 100*_ncomp*ln(1+y0/(100*_ncomp)); endif; endif; odate = getdate(tdate); if floor(tdate) < 2000; @ Para datos anteriores al año 2000 @ odate = (odate[1]-1900)~odate[2:3]; endif; if floor(tdate) >= 2000; @ Para datos desde año 2000 @ odate = (odate[1]-2000)~odate[2:3]; endif; tdstring = datestr4(odate); "*****************************************************************************"; format /rd 10,6; "COL.OUT " datestr3(0) " " timestr(0) " Trade Date: " tdstring; "*****************************************************************************"; if _yield; "Min yield errors"; else; "Min price errors"; endif; if _ncomp /=1; "Coupon and compounding frequency:" _ncomp "/yr"; endif; if restr; "Restriction: O/n rate =" _s0 "%/yr continuously compounded"; " =" 100*(exp(_s0/100)-1) "%/yr annually compounded"; if y0 == 0; " = actual o/n rate"; endif; endif; ?; "Starting parameters:"; if nsx; " b1 b2 b3 b4 b5 b6"; else; " b1 b2 b3 b4"; endif; if nsx; if restr; bstart = bst[1]|_s0-bst[1]|bst[2 3 4 5]; else; bstart = bst; endif; bstart'; ?; else; if restr; bstart = bst[1]|_s0-bst[1]|bst[2 3]; else; bstart = bst[1:4]; endif; bstart'; ?; endif; @==============================================================================@ @** Procedimientos para el manejo de Fechas **@ @==============================================================================@ /***************************************************************************** GETTIME.G 931013 LS Purpose: Constructs time measured in years from number on form yyyy.mmdd, assuming 365 days per year Format: y = GETTIME(x); Input: x NxK vector, numbers on form yyyy.mmdd Output: y NxK vector, numbers on form yyyy + (mm-1)/12 + (dd-1)/365 *****************************************************************************/ proc (1) = GETTIME(x); local yyyy,mm,dd,y; yyyy = floor(x); mm = floor( 100*(x-yyyy) ); dd = 10000*(x - yyyy - mm/100); y = yyyy + (mm-1)/12 + (dd-1)/365; retp(y); endp; /****************************************************************************/ @==============================================================================@ /***************************************************************************** GETDATE.G 930515 LS Purpose: Constructs date from number on form yyyy.mmdd Format: y = GETDATE(x); Input: x Nx1 vector, numbers on form yyyy.mmdd Output: y Nx3 vector, y[.,1] = yyyy, y[.,2] = mm, y[.,3] = dd *****************************************************************************/ proc (1) = GETDATE(x); local yyyy,mm,dd,y; yyyy = floor(x); mm = floor( 100*(x-yyyy) ); dd = 10000*(x - yyyy - mm/100); y = yyyy~mm~dd; retp(y); endp; /****************************************************************************/ @==============================================================================@ /**************************************************************************** ** DATESTR4.G ** 931207 LS ** Purpose: Formats a date in a vector to a string. ** Format: d = DATESTR3(d); ** Input: t -- 4x1 vector from the DATE function or a zero. If ** the input is 0 the DATE function will be called ** to return the current system date. ** Ouput: d -- string containing a date in the format: ** dd Mon yyyy ** Globals: None ** Note: Revision of DATESTR ** See Also: DATESTR, DATE, TIME, TIMESTR, ETHSEC, ETSTR; *****************************************************************************/ proc datestr4(dt); local m; if dt == 0; dt = DATE; endif; if dt[2] == 1; m = "Jan"; endif; if dt[2] == 2; m = "Feb"; endif; if dt[2] == 3; m = "Mar"; endif; if dt[2] == 4; m = "Apr"; endif; if dt[2] == 5; m = "May"; endif; if dt[2] == 6; m = "Jun"; endif; if dt[2] == 7; m = "Jul"; endif; if dt[2] == 8; m = "Aug"; endif; if dt[2] == 9; m = "Sep"; endif; if dt[2] == 10; m = "Oct"; endif; if dt[2] == 11; m = "Nov"; endif; if dt[2] == 12; m = "Dec"; endif; retp( ftos(dt[3],"%*.*lf",2,0) $+ " " $+ m $+ " " $+ ftos(dt[1],"%0*.*lf",2,0) ); endp; /****************************************************************************/ @==============================================================================@ /**************************************************************************** ** DATESTR3.G ** Purpose: Formats a date in a vector to a string. ** Format: d = DATESTR3(d); ** Input: t -- 4x1 vector from the DATE function or a zero. If ** the input is 0 the DATE function will be called ** to return the current system date. ** Ouput: d -- 8 character string containing a date in the format: ** yy/mm/dd ** Globals: None ** Note: Revision of DATESTR ** See Also: DATESTR, DATE, TIME, TIMESTR, ETHSEC, ETSTR; *****************************************************************************/ proc datestr3(dt); if dt == 0; dt = DATE; endif; retp( ftos(dt[1]%100,"%0*.*lf",2,0) $+ ftos(dt[2],"%0*.*lf",2,0) $+ ftos(dt[3],"%0*.*lf",2,0) ); endp; /****************************************************************************/ @==============================================================================@ /**************************************************************************** ** MISSVS.G ** LS 901220 ** Purpose: Creates a matrix of missing values ** Format: y = MISSVS(r,c); ** Input: r integer number of rows ** c integer number of columns ** Output: y rxc matrix of missing values **---------------------------------------------------------------------------*/ proc (1) = MISSVS(r,c); local y; y = ones(r,c); y = missex(y,y); retp(y); endp; /*****************************************************************************/ @==============================================================================@ /***************************************************************************** IGETTIME.G LFM 20010604 Col. Purpose: Constructs a date on form yyyy~mm~dd from time measured in years (form yyyy.mmdd) (assuming 365 days per year) Format: y = IGETTIME(x); Input: x NxK vector, numbers on form yyyy + (mm-1)/12 + (dd-1)/365 Output: y NxK vector, numbers on form yyyy~mm~dd *****************************************************************************/ proc (1) = IGETTIME(x); local yyyy,mm,dd,y; yyyy = floor(x); mm = floor( 12*(x-yyyy)+1 ); dd = int(365*(x - yyyy - (mm-1)/12)+1); y = yyyy*10000+mm*100+dd; retp(y); endp; /****************************************************************************/ @==============================================================================@ @==============================================================================@ @** Partes del programa NS2 utilizado para ML **@ @==============================================================================@ /**************************************************************************** CFTPM Purpose: Create cash flow and time to payment matrix Format: {cf,tp} = CFTPM(c,m); Input: c Nx1 vector coupons, %/yr m Nx1 vector time to maturity, yrs Output: cf NxK matrix cash flows, padded to the right with zeros tp NxK matrix times to payment, padded to the right with 99.99 cf[i,j] is cash flow nr j of bond nr i cf[i,j] = c/_ncomp if cf[i,j+1] > 0 cf[i,j] = 100 + c/_ncomp if cf[i,j+1] = 0 tp[i,j] is time to payment nr j of bond nr i Global: _ncomp scalar coupon frequency, /yr Remark: ----------------------------------------------------------------------------*/ proc (2) = CFTPM(c,m); local n,h,np,k,cf,tp,j; if rows(c) /= rows(m); "Error in c or m"; end; endif; n = rows(m); h = _ncomp; np = floor(m*h) + (m*h .> floor(m*h)); k = maxc(np); cf = zeros(n,k); tp = 99.99*ones(n,k); j = 1; do while j <= n; if np[j] > 1; cf[j,1:(np[j]-1)] = (c[j]/h)*ones(1,np[j]-1); endif; cf[j,np[j]] = 100 + c[j]/h; tp[j,1:np[j]] = m[j] - floor(m[j]*h)/h + ( 1/h - (m[j] > floor(m[j]*h)/h)/h ) + seqa(0,1/h,np[j])'; j = j+1; endo; retp(cf,tp); endp; @==============================================================================@ /************************************************************************** PV2 Purpose: Generate present values of coupon bonds Format: p = PV2(y); Input: y Nx1 vector yields to maturity, %/yr Output: p Nx1 vector prices of coupon bonds Globals: _ncomp scalar coupon and compounding frequency, /yr _cf NxK matrix cash flows (from CFTPM) _tp NxK matrix times to payment (from CFTPM) Remark: ----------------------------------------------------------------------------*/ proc PV2(y); local h,p; h = _ncomp; p = diag( _cf * (1 + y'/(100*h))^(-h*_tp') ); retp(p); endp; @==============================================================================@ /**************************************************************************** SPOTR Purpose: Nelson and Siegel (1987) spot rate function. Format: s = SPOTR(m,b); Input: m NxK matrix times to maturity, yrs b 4x1 vector b[1] = 100*beta0 in NS, %/yr b[2] = 100*beta1 in NS, %/yr b[3] = 100*beta2 in NS, %/yr b[4] = tau in NS, yr Output: s NxK matrix spot rate, %/yr Remark: Spot rate continously compounded at annual rate. ----------------------------------------------------------------------------*/ proc SPOTR(m,b); local s; s = b[1]/100 + (b[2]/100+b[3]/100)*(1 - exp(-m/b[4]))./(m/b[4]) - (b[3]/100)*exp(-m/b[4]); retp( 100*s ); endp; @==============================================================================@ /**************************************************************************** FORWR Purpose: Nelson and Siegel (1987) forward rate function. Format: f = FORWR(m,b); Input: m NxK matrix times to maturity, yrs b 4x1 vector b[1] = 100*beta0 in NS, %/yr b[2] = 100*beta1 in NS, %/yr b[3] = 100*beta2 in NS, %/yr b[4] = tau in NS, yr Output: f NxK matrix instantaneous forward rate, %/yr Remark: Forward rate continouosly compounded at annual rate. ----------------------------------------------------------------------------*/ proc FORWR(m,b); local f; f = b[1]/100 + (b[2]/100)*exp(-m/b[4]) + (b[3]/100)*(m/b[4]).*exp(-m/b[4]); retp( 100*f ); endp; @==============================================================================@ /**************************************************************************** DISCFN Purpose: Nelson and Siegel (1987) discount function. Format: d = DISCFN(m,b); Input: m NxK matrix times to maturity, yrs b 4x1 vector b[1] = 100*beta0 in NS, %/yr b[2] = 100*beta1 in NS, %/yr b[3] = 100*beta2 in NS, %/yr b[4] = tau in NS, yr Output: d NxK matrix discount factors, prices of unit discount bonds. Remark: Uses proc SPOTR. ----------------------------------------------------------------------------*/ proc DISCFN(m,b); external proc SPOTR; local s,d; s = SPOTR(m,b); d = exp( -(s/100).*m ); retp( d ); endp; @==============================================================================@ /*************************************************************************** PVNS2 Purpose: Generate NS present values of coupon bonds Format: p = PVNS2(b); Input: b 4x1 vector parameters b[1] is 100*beta0 in NS (1987), %/yr b[2] is 100*beta1 in NS, %/yr b[3] is 100*beta2 in NS, %/yr b[4] is tau in NS, yr Output: p Nx1 vector prices of coupon bonds Globals: _cf NxK matrix cash flows (from CFTPM) _tp NxK matrix times to payment (from CFTPM) Remark: ----------------------------------------------------------------------------*/ proc PVNS2(b); external proc DISCFN; local p; b[4] = abs(b[4]); /* Avoid crash if negative tau */ p = diag(_cf*DISCFN(_tp',b)); retp(p); endp; @==============================================================================@ /************************************************************************** LLIK2 Purpose: Log likelihood funktion Format: ll = LLIK2(b,pdat); Input: b 4x1 vector parameters b[1] is 100*beta0 in NS (1987), %/yr b[2] is 100*beta1 in NS, %/yr b[2] is 100*beta2 in NS, %/yr b[3] is tau in NS, yr pdat Nx1 matrix bond price data Output: ll Nx1 vector log likelihoods Globals: _cf NxK matrix cash flows (from CFTPM) _tp NxK matrix times to payment (from CFTPM) Remark: Prevents crash for negative taus. ----------------------------------------------------------------------------*/ proc LLIK2(b,pdat); local bb,dev,s2; bb = b; bb[4] = abs(bb[4]); /* Prevents crash */ dev = (PVNS2(bb) - pdat); s2 = dev'dev/rows(dev); retp(-0.5*(dev.*dev/s2 + ln(2*pi*s2))); endp; @==============================================================================@ /************************************************************************** LLIK2X Purpose: Log likelihood funktion Format: ll = LLIK2X(b,pdat); Input: b 6x1 vector b[1] = 100*beta0 in NS, %/yr b[2] = 100*beta1 in NS, %/yr b[3] = 100*beta2 in NS, %/yr b[4] = tau in NS, yr b[5] = 100*beta3, %/yr b[6] = tau2, yr pdat Nx1 matrix bond price data Output: ll Nx1 vector log likelihoods Globals: _cf NxK matrix cash flows (from CFTPM) _tp NxK matrix times to payment (from CFTPM) Remark: Avoids crash for negative taus. ----------------------------------------------------------------------------*/ proc LLIK2X(b,pdat); local bb,dev,s2; bb = b; bb[4] = abs(bb[4]); /* Prevent crash for negative taus */ bb[6] = abs(bb[6]); dev = (PVNS2X(bb) - pdat); s2 = dev'dev/rows(dev); retp(-0.5*(dev.*dev/s2 + ln(2*pi*s2))); endp; @==============================================================================@ /*************************************************************************** LLIK21 Purpose: Log likelihood funktion, with given o/n rate Format: ll = LLIK21(b,pdat); Input: b 3x1 vector parameters b[1] is 100*beta0 in NS (1987), %/yr b[2] is 100*beta2 in NS, %/yr b[3] is tau in NS, yr _s0-b[1] is 100*beta1 in NS, %/yr pdat Nx1 vector bond price data Output: ll Nx1 vector log likelihoods Globals: _cf NxK matrix cash flows (from CFTPM) _tp NxK matrix times to payment (from CFTPM) _s0 scalar o/n rate, %/yr Remark: ----------------------------------------------------------------------------*/ proc LLIK21(b,pdat); local bb,dev,s2; bb = b[1]|(_s0-b[1])|b[2]|b[3]; bb[4] = abs(bb[4]); /* Prevents crash */ dev = (PVNS2(bb) - pdat); s2 = dev'dev/rows(dev); retp(-0.5*(dev.*dev/s2 + ln(2*pi*s2))); endp; @==============================================================================@ /*************************************************************************** LLIK21X Purpose: Log likelihood funktion, with given o/n rate Format: ll = LLIK21X(b,pdat); Input: b 5x1 vector parameters b[1] is 100*beta0 in NS (1987), %/yr b[2] is 100*beta2 in NS, %/yr b[3] is tau in NS, yr b[4] = 100*beta3, %/yr b[5] = tau2, yr _s0-b[1] is 100*beta1 in NS, %/yr pdat Nx1 vector bond price data Output: ll Nx1 vector log likelihoods Globals: _cf NxK matrix cash flows (from CFTPM) _tp NxK matrix times to payment (from CFTPM) _s0 scalar o/n rate, %/yr Remark: ----------------------------------------------------------------------------*/ proc LLIK21X(b,pdat); local bb,dev,s2; bb = b[1]|(_s0-b[1])|b[2]|b[3]|b[4]|b[5]; bb[4] = abs(bb[4]); /* Prevent crash for negative taus */ bb[6] = abs(bb[6]); dev = (PVNS2X(bb) - pdat); s2 = dev'dev/rows(dev); retp(-0.5*(dev.*dev/s2 + ln(2*pi*s2))); endp; @==============================================================================@ /**************************************************************************** YTM2 Purpose: Compute yields to maturity from prices, with Newton-Raphson algorithm Format: ytm = YTM2(p); Input: p Nx1 vector prices Output: ytm Nx1 vector annualized yields to maturity, compounded _ncomp per year, %/yr Globals: _ncomp scalar compounding frequency _cf NxK matrix cash flows (from CFTPM) _tp NxK matrix times to payment (from CFTPM) _ytol scalar tolerance Remark: Input is column vector, not row vectors as in YAC ----------------------------------------------------------------------------*/ proc YTM2(p); local h,Dp,y00,y0,eps,y; h = _ncomp; y00 = 100*h*( (_cf[.,1]./p)^(1./(_tp[.,1]*h)) - 1 ); y0 = (_cf[.,2] .== 0).*y00 + (_cf[.,2] .> 0).*_cf[.,1]*h; eps = 1; do while eps > _ytol; _ptest = p; /* For testing */ _ytest = y0; /* if minc(y0) < 0; end; endif; */ /* "_ptest = " _ptest'; " _ytest = " _ytest'; */ Dp = diag(_cf * ((-_tp'*h).*(1+y0'/(100*h))^(-_tp'*h-1)) /(100*h)); y = y0 + (_cf[.,2] .> 0).*(p-PV2(y0))./Dp; eps = maxc( abs(y-y0) ); /* "eps =" eps; */ y0 = y; y0 = abs(y); /* to prevent crash */ endo; retp( y ); endp; @==============================================================================@ /************************************************************************** LLIK3 Purpose: Log likelihood funktion Format: ll = LLIK3(b,ydat); Input: b 4x1 vector parameters b[1] is 100*beta0 in NS (1987), %/yr b[2] is 100*beta1 in NS, %/yr b[2] is 100*beta2 in NS, %/yr b[3] is tau in NS, yr ydat Nx1 matrix yield to maturity data Output: ll Nx1 vector log likelihoods Globals: _cf NxK matrix cash flows (from CFTPM) _tp NxK matrix times to payment (from CFTPM) _ytol scalar tolerance for YTM Remark: ----------------------------------------------------------------------------*/ proc LLIK3(b,ydat); local bb,dev,s2; bb = b; bb[4] = abs(bb[4]); /* Prevents crash */ dev = (YTM2(PVNS2(bb)) - ydat); s2 = dev'dev/rows(dev); retp(-0.5*(dev.*dev/s2 + ln(2*pi*s2))); endp; @==============================================================================@ /************************************************************************** LLIK3X Purpose: Log likelihood funktion Format: ll = LLIK3X(b,ydat); Input: b 6x1 vector parameters b[1] is 100*beta0 in NS (1987), %/yr b[2] is 100*beta1 in NS, %/yr b[3] is 100*beta2 in NS, %/yr b[4] is tau in NS, yr b[5] = 100*beta3, %/yr b[6] = tau2, yr ydat Nx1 matrix yield to maturity data Output: ll Nx1 vector log likelihoods Globals: _cf NxK matrix cash flows (from CFTPM) _tp NxK matrix times to payment (from CFTPM) _ytol scalar tolerance for YTM Remark: ----------------------------------------------------------------------------*/ proc LLIK3X(b,ydat); local bb,dev,s2; bb = b; bb[4] = abs(bb[4]); /* Prevent crash for negative taus */ bb[6] = abs(bb[6]); dev = (YTM2(PVNS2X(bb)) - ydat); s2 = dev'dev/rows(dev); retp(-0.5*(dev.*dev/s2 + ln(2*pi*s2))); endp; @==============================================================================@ /*************************************************************************** LLIK31 Purpose: Log likelihood funktion, with given o/n rate Format: ll = LLIK31(b,ydat); Input: b 3x1 vector parameters b[1] is 100*beta0 in NS (1987), %/yr b[2] is 100*beta2 in NS, %/yr b[3] is tau in NS, yr _s0-b[1] is 100*beta1 in NS, %/yr ydat Nx1 vector yield to maturity data Output: ll Nx1 vector log likelihoods Globals: _cf NxK matrix cash flows (from CFTPM) _tp NxK matrix times to payment (from CFTPM) _s0 scalar o/n rate, %/yr _ytol scalar tolerance for YTM Remark: ----------------------------------------------------------------------------*/ proc LLIK31(b,ydat); local bb,dev,s2; bb = b[1]|(_s0-b[1])|b[2]|b[3]; bb[4] = abs(bb[4]); /* Prevents crash */ dev = (YTM2(PVNS2(bb)) - ydat); s2 = dev'dev/rows(dev); retp(-0.5*(dev.*dev/s2 + ln(2*pi*s2))); endp; @==============================================================================@ /*************************************************************************** LLIK31X Purpose: Log likelihood funktion, with given o/n rate Format: ll = LLIK31X(b,ydat); Input: b 5x1 vector parameters b[1] is 100*beta0 in NS (1987), %/yr b[2] is 100*beta2 in NS, %/yr b[3] is tau in NS, yr b[4] = 100*beta3, %/yr b[5] = tau2, yr _s0-b[1] is 100*beta1 in NS, %/yr ydat Nx1 vector yield to maturity data Output: ll Nx1 vector log likelihoods Globals: _cf NxK matrix cash flows (from CFTPM) _tp NxK matrix times to payment (from CFTPM) _s0 scalar o/n rate, %/yr _ytol scalar tolerance for YTM Remark: ----------------------------------------------------------------------------*/ proc LLIK31X(b,ydat); local bb,dev,s2; bb = b[1]|(_s0-b[1])|b[2]|b[3]|b[4]|b[5]; bb[4] = abs(bb[4]); /* Prevent crash for negative taus */ bb[6] = abs(bb[6]); dev = (YTM2(PVNS2X(bb)) - ydat); s2 = dev'dev/rows(dev); retp(-0.5*(dev.*dev/s2 + ln(2*pi*s2))); endp; @==============================================================================@ /*------------------------------------------------------------------------------ NSEST2 Purpose: Nelson and Siegel (1987) estimate of term structure Format: { b,fmin,g,cov,retcode,rmspe,rmsye } = NSEST2(c,m,y,bst); Input: c Nx1 vector coupons, %/yr double semiannual coupons if semiannual paym. m Nx1 vector times to maturity, years y Nx1 vector yields to maturity, %/yr bst 4x1 vector starting values, bst[1:3] %/yr Output: b 4x1 matrix estiamated parameters, b[1:3] %/yr b[1] is 100*beta0 in NS, %/yr b[2] is 100*beta1 in NS, %/yr b[3] is 100*beta2 in NS, %/yr b[4] is tau in NS, yr cov 4x4 matrix estimated covariance of b if _s0 = 0 3x3 matrix estimated cov of b[1 3 4] if _s0r > 0 fmin scalar minimum log likelihood g vector gradient at min retcode scalar return code from MAXLIK rmspe scalar root mean squared price error, % of face value rmsye scalar root mean squared yield error, %/yr _cf NxK matrix cash flows _tp NxK matrix times to payment Globals: _ncomp scalar compounding and coupon frequence, /yr _s0 scalar ==0 unrestricted estimate >0 given o/n rate, %/yr _yield scalar 0= min price error, 1= min yield error _ytol scalar tolerance for YTM _cf cash flows, _cf = 0 to initilize _tp times to payment, _tp = 0 to initialize Remark: If _s0 > 0 only bst[1 3 4] used. ------------------------------------------------------------------------------*/ proc (7) = NSEST2(c,m,y,bst); external matrix _ncomp,_s0,_cf,_tp; local pdat,bb,fmin,g,cov,retcode,phat,rmspe,yhat,rmsye,b,bst1,b1; {_cf,_tp} = CFTPM(c,m); pdat = PV2(y); /*------------------------------------------------------------------------- ML Optimization --------------------------------------------------------------------------*/ maxset; if max_algorit /= 0; _max_Algorithm = max_algorit; @ 2->BFGLS, 3->DFP, 4->NEWTON; @ endif; if max_mlcov /= 0; _max_CovPar = max_mlcov; @ 0->Not computed, 1->from Hessian(def) , 2->from QML; @ endif; if diagnost_ml == 1; _max_Diagnostic = 1; endif; _max_GradTol = toler; @_max_GradTol = 1e-5;@ @Default;@ __output = 0; if not _yield; if _s0 == 0; __title = "Min price error, w/o restriction "; {bb,fmin,g,cov,retcode} = maxlik(pdat,1,&llik2,bst); else; __title = "Min price error, w/ restriction "; bst1 = bst[1 3 4]; {bb,fmin,g,cov,retcode} = maxlik(pdat,1,&llik21,bst1); endif; else; if _s0 == 0; __title = "Min yield error, w/o restriction "; {bb,fmin,g,cov,retcode} = maxlik(y,1,&llik3,bst); else; __title = "Min yield error, w/ restriction "; bst1 = bst[1 3 4]; {bb,fmin,g,cov,retcode} = maxlik(y,1,&llik31,bst1); endif; endif; call maxprt(bb,fmin,g,cov,retcode); if _s0 == 0; b = bb; else; b = bb[1]|_s0-bb[1]|bb[2 3]; endif; b[4] = abs(b[4]); /* Prevents crash */ phat = PVNS2(b); rmspe = sqrt((phat-pdat)'(phat-pdat)/rows(pdat)); yhat = YTM2(phat); rmsye = sqrt((yhat-y)'(yhat-y)/rows(pdat)); retp( b,fmin,g,cov,retcode,rmspe,rmsye); endp; @==============================================================================@ /**************************************************************************** NSEST2X Purpose: Nelson and Siegel (1987) estimate of term structure Format: { b,fmin,g,cov,retcode,rmspe,rmsye } = NSEST2X(c,m,y,bst); Input: c Nx1 vector coupons, %/yr double semiannual coupons if semiannual paym. m Nx1 vector times to maturity, years y Nx1 vector yields to maturity, %/yr bst 4x1 vector starting values, bst[1:3] %/yr Output: b 6x1 matrix estiamated parameters, b[1:3] %/yr b[1] is 100*beta0 in NS, %/yr b[2] is 100*beta1 in NS, %/yr b[3] is 100*beta2 in NS, %/yr b[4] is tau in NS, yr b[5] is 100*beta2, %/yr b[6] is tau2, yr cov 6x6 matrix estimated covariance of b if _s0 = 0 5x5 matrix estimated cov of b[1 3 4 5 6] if _s0r > 0 fmin scalar minimum log likelihood g vector gradient at min retcode scalar return code from MAXLIK rmspe scalar root mean squared price error, % of face value rmsye scalar root mean squared yield error, %/yr _cf NxK matrix cash flows _tp NxK matrix times to payment Globals: _ncomp scalar compounding and coupon frequence, /yr _s0 scalar ==0 unrestricted estimate >0 given o/n rate, %/yr _yield scalar 0= min price error, 1= min yield error _ytol scalar tolerance for YTM _cf cash flows, _cf = 0 to initilize _tp times to payment, _tp = 0 to initialize Remark: If _s0 > 0 only bst[1 3 4 5 6] used. ------------------------------------------------------------------------------*/ proc (7) = NSEST2X(c,m,y,bst); external matrix _ncomp,_s0,_cf,_tp; local pdat,bb,fmin,g,cov,retcode,phat,rmspe,yhat,rmsye,b, bst1,b1; {_cf,_tp} = CFTPM(c,m); pdat = PV2(y); /*------------------------------------------------------------------------- ML Optimization --------------------------------------------------------------------------*/ maxset; /* _mlcovp = 3;*/ @ __max_Algorithm = 2 ==> BFGLS , Default: 2 ==> DFP @ _max_Algorithm = 2; _max_CovPar = 2; __output = 0; if not _yield; if _s0 == 0; __title = "Min price error, w/o restriction "; {bb,fmin,g,cov,retcode} = maxlik(pdat,1,&llik2x,bst); else; __title = "Min price error, w/ restriction "; bst1 = bst[1 3 4 5 6]; {bb,fmin,g,cov,retcode} = maxlik(pdat,1,&llik21x,bst1); endif; else; if _s0 == 0; __title = "Min yield error, w/o restriction "; {bb,fmin,g,cov,retcode} = maxlik(y,1,&llik3x,bst); else; __title = "Min yield error, w/ restriction "; bst1 = bst[1 3 4 5 6]; {bb,fmin,g,cov,retcode} = maxlik(y,1,&llik31x,bst1); endif; endif; call maxprt(bb,fmin,g,cov,retcode); if _s0 == 0; b = bb; else; b = bb[1]|_s0-bb[1]|bb[2 3 4 5]; endif; b[4] = abs(b[4]); b[6] = abs(b[6]); /* Correct negative taus */ phat = PVNS2X(b); rmspe = sqrt((phat-pdat)'(phat-pdat)/rows(pdat)); yhat = YTM2(phat); rmsye = sqrt((yhat-y)'(yhat-y)/rows(pdat)); retp( b,fmin,g,cov,retcode,rmspe,rmsye); endp; /******************************************************************************/ @==============================================================================@ /******************************************************************************* ML Optimization *******************************************************************************/ if opt; _cf = 0; _tp = 0; /* These two globals must be initialized */ if nsx; {bh,fmin,g,cov,retcode,rmspe,rmsye} = NseST2X(c',tm',y',bstart); else; {bh,fmin,g,cov,retcode,rmspe,rmsye} = NseST2(c',tm',y',bstart); bh = bh | 0 | 1; endif; endif; /*****************************************************************************/ pdat = PV2(y'); @==============================================================================@ /****************************************************************************** SPOTRX Purpose: Extended Nelson and Siegel (1987) spot rate function. Format: s = SPOTRX(m,b); Input: m NxK matrix times to maturity, yrs b 6x1 vector b[1] = 100*beta0 in NS, %/yr b[2] = 100*beta1 in NS, %/yr b[3] = 100*beta2 in NS, %/yr b[4] = tau in NS, yr b[5] = 100*beta3, %/yr b[6] = tau2, yr Output: s NxK matrix spot rate, %/yr Remark: Spot rate continously compounded at annual rate. ----------------------------------------------------------------------------*/ proc SPOTRX(m,b); local s; s = b[1]/100 + (b[2]/100)*(1 - exp(-m/b[4]))./(m/b[4]) + (b[3]/100)*((1 - exp(-m/b[4]))./(m/b[4]) - exp(-m/b[4])) + (b[5]/100)*((1 - exp(-m/b[6]))./(m/b[6]) - exp(-m/b[6])); retp( 100*s ); endp; @==============================================================================@ /**************************************************************************** FORWRX Purpose: Extended Nelson and Siegel (1987) forward rate function. Format: f = FORWRX(m,b); Input: m NxK matrix times to maturity, yrs b 6x1 vector b[1] = 100*beta0 in NS, %/yr b[2] = 100*beta1 in NS, %/yr b[3] = 100*beta2 in NS, %/yr b[4] = tau in NS, yr b[5] = 100*beta3, %/yr b[6] = tau2, yr Output: f NxK matrix instantaneous forward rate, %/yr Remark: Forward rate continouosly compounded at annual rate. ----------------------------------------------------------------------------*/ proc FORWRX(m,b); local f; f = b[1]/100 + (b[2]/100)*exp(-m/b[4]) + (b[3]/100)*(m/b[4]).*exp(-m/b[4]) + (b[5]/100)*(m/b[6]).*exp(-m/b[6]); retp( 100*f ); endp; @==============================================================================@ /**************************************************************************** DISCFNX Purpose: Extended Nelson and Siegel (1987) discount function. Format: d = DISCFNX(m,b); Input: m NxK matrix times to maturity, yrs b 6x1 vector b[1] = 100*beta0 in NS, %/yr b[2] = 100*beta1 in NS, %/yr b[3] = 100*beta2 in NS, %/yr b[4] = tau in NS, yr b[5] = 100*beta3, %/yr b[6] = tau2, yr Output: d NxK matrix discount factors, prices of unit discount bonds. Remark: Uses proc SPOTRX. ----------------------------------------------------------------------------*/ proc DISCFNX(m,b); external proc SPOTRX; local s,d; s = SPOTRX(m,b); d = exp( -(s/100).*m ); retp( d ); endp; @==============================================================================@ /*************************************************************************** PVNS2X Purpose: Generate extended NS present values of coupon bonds Format: p = PVNS2X(b); Input: b 6x1 vector b[1] = 100*beta0 in NS, %/yr b[2] = 100*beta1 in NS, %/yr b[3] = 100*beta2 in NS, %/yr b[4] = tau in NS, yr b[5] = 100*beta3, %/yr b[6] = tau2, yr Output: p Nx1 vector prices of coupon bonds Globals: _cf NxK matrix cash flows (from CFTPM) _tp NxK matrix times to payment (from CFTPM) Remark: ----------------------------------------------------------------------------*/ proc PVNS2X(b); external proc DISCFNX; local p; b[4] = abs(b[4]); /* Avoid crash for negative tau and tau2 */ b[6] = abs(b[6]); p = diag(_cf*DISCFNX(_tp',b)); retp(p); endp; @==============================================================================@ /******************************************************************************/ /* Estimates */ ph = PVNS2X(bh); /* Coupon bond prices */ sh = 100*(exp(SPOTRX(tm',bh)/100) - 1);/* Annually compounded spot rate */ sinfh = 100*(exp(bh[1]/100) - 1);/* Annually compounded asympt spot rate */ fh = 100*(exp(FORWRX(tm',bh)/100) - 1);/*Annually compounded forward rate*/ yh = YTM2(ph); /******************************************************************************/ /* Plot curves */ @ tmstep = (tm[nm] - mmin)/(nm1-1); @ temp13 = sortc(tm',1); tmstep = (temp13[nm] - mmin)/(nm1-1); tm1 = seqa(mmin,tmstep,nm1)'; ph1 = 100*DISCFNX(tm1,bh); /* Discount bond prices */ sh1 = 100*(exp(SPOTRX(tm1,bh)/100) - 1);/* Annually compounded spot rate */ fh1 = 100*(exp(FORWRX(tm1,bh)/100)-1); /*Annually comp forward rate */ tmp = tm~missvs(1,nm1-nm); /* Pad for plot */ yp = y~missvs(1,nm1-nm); /* Pad for plot */ yhp = yh'~missvs(1,nm1-nm); /* Pad for plot */ cp = c~missvs(1,nm1-nm); /* Pad for plot */ /******************************************************************************/ /* Error bars */ tm2 = seqa(ebmin,ebstep,neb); tm2 = tm2 + ebfact1*tm2^ebfact2; proc SPOTR2X(b2); local bb,s; bb = b2[1]|_s0-b2[1]|b2[2 3 4 5]; s = SPOTRX(tm2,bb); retp(s); endp; proc FORWR2X(b2); local bb,f; bb = b2[1]|_s0-b2[1]|b2[2 3 4 5]; f = FORWRX(tm2,bb); retp( f ); endp; proc PVNS22X(b2); local bb,p; bb = b2[1]|_s0-b2[1]|b2[2 3 4 5]; p = PVNS2X(bb); retp( p ); endp; proc YTM22X(b2); local bb,pp,yy; bb = b2[1]|_s0-b2[1]|b2[2 3 4 5]; pp = PVNS2X(bb); yy = YTM2(pp); retp( yy ); endp; proc SPOTR3X(b); local s; s = SPOTRX(tm2,b); retp(s); endp; proc FORWR3X(b); local f; f = FORWRX(tm2,b); retp( f ); endp; proc YTM3X(b); retp( YTM2(PVNS2X(b)) ); endp; proc SPOTR2(b2); local bb,s; bb = b2[1]|_s0-b2[1]|b2[2 3]; s = SPOTR(tm2,bb); retp(s); endp; proc FORWR2(b2); local bb,f; bb = b2[1]|_s0-b2[1]|b2[2 3]; f = FORWR(tm2,bb); retp( f ); endp; proc PVNS22(b2); local bb,p; bb = b2[1]|_s0-b2[1]|b2[2 3]; p = PVNS2(bb); retp( p ); endp; proc YTM22(b2); local bb,pp,yy; bb = b2[1]|_s0-b2[1]|b2[2 3]; pp = PVNS2(bb); yy = YTM2(pp); retp( yy ); endp; proc SPOTR3(b); local s; s = SPOTR(tm2,b); retp(s); endp; proc FORWR3(b); local f; f = FORWR(tm2,b); retp( f ); endp; proc YTM3(b); retp( YTM2(PVNS2(b)) ); endp; /******************************************************************************/ format /rd 10,6;?; "Optimal parameters:"; if nsx; " b1 b2 b3 b4 b5 b6"; else; " b1 b2 b3 b4"; endif; if nsx; bh'; ?; else; bh[1:4]'; ?; endif; "Infinite maturity spotrate:" sinfh "%/yr annually compounded";?; "Price errors:"; " Mean squared price error:" rmspe^2 ; " RMSPE:" rmspe; " Mean absolute price error:" meanc(abs(ph-pdat)); " Max absolute price error:" maxc(abs(ph-pdat)); tmm = tm[maxindc(abs(ph-pdat))]; " for time to maturity:" tmm "yrs = " 12*tmm "months";?; msye = (yh-y')'(yh-y') / nm; "Yield errors"; " Mean squared yield error:" msye "%/yr"; " RMSYE:" sqrt(msye) "%/yr"; " Mean absolute yield error:" meanc(abs(yh-y')) "%/yr"; " Max absolute yield error:" maxc(abs(yh-y')) "%/yr"; tmm = tm[maxindc(abs(yh-y'))]; " for time to maturity:" tmm "yrs = " 12*tmm "months"; @" ";@ @" ";@ @pdat';@ @*" ";@ /******************************************************************************/ /* Compute elapsed time */ etime1 = ethsec(tstart,date); etime = etstr(etime1); "*****************************************************************************"; "Program completed after " etime; "*****************************************************************************"; @fname = "c:\\archivos\\excel\\temp.xls";@ @let names = b1 b2 b3 b4 std_b2 std_b3 std_b4 rmspe_p mape_p rmspe_y mape_y;@ @x111 = bh[1:4]'~sqrt(diag(cov[1:3,1:3]))'~rmspe~meanc(abs(ph-pdat))~sqrt(msye)~meanc(abs(yh-y'));@ @call export(x,fname,names);@ @save c:\archivos\excel\temp = x111;@ @Salida utilizada para el diagnóstico para N-S con restricciones @ @bh[1:4]'~sqrt(diag(cov[1:3,1:3]))'~rmspe~meanc(abs(ph-pdat))~sqrt(msye)~meanc(abs(yh-y'));@ @Salida utilizada para el diagnóstico para N-S sin O/N sin restricciones @ @bh[1:4]'~sqrt(diag(cov[1:4,1:4]))'~rmspe~meanc(abs(ph-pdat))~sqrt(msye)~meanc(abs(yh-y'));@ output off; if prf; dos lp c:\archivos\term_str\svensson\col\COLX.out /L5 /W95 /H85; endif; @==============================================================================@ /***************************************************************************** CONFNF.G 931003 Purpose: Computes normal confidence interval of function of estimates by first-order Taylor expansion, the delta method. Format: y = CONFNF(&f,b,cov,q,tol); Input: f proc mapping b into Nx1 vector b Kx1 vector estimates cov KxK matrix variance matrix of estimates q scalar confidence level tol scalar tolerance Output: y Nx3 vector y[.,1] lower limit of confidence interval y[.,2] point estimate of linear combination y[.,3] upper limit of confidence interval Remark: Uses external proc CDFNINV. *****************************************************************************/ proc (1) = CONFNF(&f,b,cov,q,tol); external proc CDFNINV; local nn,yy,fb,f:proc,sey,y; nn = CDFNINV((1+q)/2,tol); yy = f(b); fb = gradp(&f,b); sey = sqrt(diag(fb*cov*fb')); y = zeros(rows(yy),3); y[.,2] = yy; y[.,1] = y[.,2] - nn*sey; y[.,3] = y[.,2] + nn*sey; retp(y); endp; /****************************************************************************/ @==============================================================================@ /***************************************************************************** CDFNINV.G 911031 Purpose: Computes inverse of normal cumulative distribution function Format: y = CDFNINV(x,tol); Input: x scalar value of cumulative distribution function tol scalar tolerance Output: y scalar inverse of CDFN, x = CDFN(y) *****************************************************************************/ proc (1) = CDFNINV(x,tol); local ymin,ymax,y0,y; if x <= 0 or x >= 1; errorlog "ERROR: First argument must be strictly between 0 and 1"; end; endif; if tol <= 0; errorlog "ERROR: tolerance must be positive"; end; endif; ymin = 0; do until cdfn(ymin) < x; ymin = ymin - 1; endo; ymax = 0; do until cdfn(ymax) > x; ymax = ymax + 1; endo; y = (ymin + ymax)/2; y0 = y + 1; do until abs(y-y0) < tol; if cdfn(y) < x; y0 = y; ymin = y; y = (y + ymax)/2; elseif cdfn(y) > x; y0 = y; ymax = y; y = (y + ymin)/2; else; break; endif; endo; retp(y); endp; @==============================================================================@ /* Plot */ graphset; /* graphprt("-w=2"); */ _pdate = ""; if prgr; graphprt("-P -PR=1"); endif; if pdate; @_pdate = "LS/MB ";@ _pdate = " "; endif; if mess; _pmsgctl = {0 6.7 .10 0 2 15 0}; @ _pmsgstr = "LS/MB " $+ datestr3(0);@ _pmsgstr = datestr3(0); endif; _pgrid = {0,0}; /* No grid */ if grid; _pgrid = {1,1}; endif; /* Grid */ _pframe = {0,0}; _plctrl = {-1,0,1,1,1}; _pnum = 2; _plwidth = 3; if plgr1; _ptek="c:\\archivos\\term_str\\svensson\\col\\colfig1.tkf"; @title(" "\ "\LYields to Maturity, Spot and Forward Rates, Colombia, " $+ tdstring);@ title(" "\ "\LTasas spot, forward y de rendimiento al vencimiento, " $+ tdstring); if errbar; @title(" "\ "\LYields to Maturity, Spot and Forward Rates, Colombia, " $+ tdstring $+ "\L95% confidence interval" );@ title(" "\ "\LTasas spot, forward y de rendimiento al vencimiento, " $+ tdstring $+ "\LIntervalos de confianza del 95%" ); endif; if restr==0 and nsx == 0; title(" "\ "\LTasas spot, forward y de rendimiento al vencimiento, " $+ tdstring $+ "\LMetodologia de Nelson y Siegel sin restriccion" ); endif; if restr==1 and nsx == 0; title(" "\ "\LTasas spot, forward y de rendimiento al vencimiento, " $+ tdstring $+ "\LMetodologia de Nelson y Siegel con restriccion" ); endif; @xlabel("Maturity/settlement");@ xlabel("Vencimiento/contrato"); @ylabel("%/yr, annually compounded");@ ylabel("Tasas compuestas anuales (%)"); if asy; @ _pline = 1~1~ttm+5~sinfh~ttm+40~sinfh~1~14~0;@ _pline = 1~1~ttm+4~sinfh~ttm+40~sinfh~1~2~0; endif; _plctrl = { -1, 0, 0,-1,-1 }; @_pcolor = { 10,11,12,13,4 };@ @_pcolor = { 10,11,12,13,13 };@ _pcolor = { 3,4,12,13,1 }; _pltype = { 4, 3, 6, 1,5 }; _pstype = { 2, 3, 1, 8,4 }; _psymsiz = {3.5, 1, 1, 1,3.5}; @_plegctl = {2,4,5.4,3.95};@ _plegctl = {2,4,4.5,3.95}; @_plegstr = "Observed yield to maturity"\ "\000Estimated spot rate"\ "\000Estimated forward rate"\ "\000Estimated yield to maturity"\ "\000Coupon rate";@ _plegstr = "Rendimiento al vencimiento observado"\ "\000Tasa spot estimada"\ "\000Tasa forward estimada"\ "\000Rendimiento al vencimiento estimado"\ "\000Tasa cupon"; if errbar; if nsx; if restr; b2 = bh[1 3 4 5 6]; ys = 100*(exp(CONFNF(&spotr2x,b2,cov,q,tol)/100)-1); yf = 100*(exp(CONFNF(&forwr2x,b2,cov,q,tol)/100)-1); yy = CONFNF(&YTM22x,b2,cov,q,tol); else; ys = 100*(exp(CONFNF(&spotr3x,bh,cov,q,tol)/100)-1); yf = 100*(exp(CONFNF(&forwr3x,bh,cov,q,tol)/100)-1); yy = CONFNF(&YTM3x,bh,cov,q,tol); endif; else; if restr; b2 = bh[1 3 4]; ys = 100*(exp(CONFNF(&spotr2,b2,cov,q,tol)/100)-1); yf = 100*(exp(CONFNF(&forwr2,b2,cov,q,tol)/100)-1); yy = CONFNF(&YTM22,b2,cov,q,tol); else; b2 = bh[1:4]; ys = 100*(exp(CONFNF(&spotr3,b2,cov,q,tol)/100)-1); yf = 100*(exp(CONFNF(&forwr3,b2,cov,q,tol)/100)-1); yy = CONFNF(&YTM3,b2,cov,q,tol); endif; endif; _perrbar = ttm+(tm2~tm2~tm2)~ys[.,2 1 3]~ones(rows(tm2),1)*(6~0~0) | ttm+(tm2~tm2~tm2)~yf[.,2 1 3]~ones(rows(tm2),1)*(6~0~0) | ttm+(tm'~tm'~tm')~yy[.,2 1 3]~ones(rows(tm'),1)*(6~0~0); endif; scale(ttm+(tmp'~tm1'~tm1'), (yp'~sh1'~fh1')); if xtics1 == 1; xtics(2000,2009,1,2); endif; if xtics1 == 2; xtics(2000,2009,1,2); endif; if ytics1 == 1; @ ytics(3,13,2,2); @ @ytics(10,23,2,2); @ ytics(9,26,2,2); endif; _pmcolor={"","","","","","","","",15}; @graphprt("-P -C=4 -PF=graf1.bit");@ graphprt("-C=4"); xy(ttm+(tmp'~tm1'~tm1'~tmp'~tmp'), (yp'~sh1'~fh1'~yhp'~cp')); endif; if output_file; @for i (1,cols(tm1),1);@ @ temp1[i] = datestr3((ttm+tm1[i]));@ @endfor;@ @temp1 = igettime(ttm+tm1');@ /* Spot, Forward and Discount function */ /* fname = "c:\\archivos\\term_str\\svensson\\col\\ts_sfd.xls"; let names = maturity espot eforward ediscount; x = (ttm+tm1')~sh1'~fh1'~ph1'; call export(x,fname,names); */ /* Yield */ /* fname = "c:\\archivos\\term_str\\svensson\\col\\ts_y.xls"; let names = maturity oyield eyield; x = (ttm+tmp')~yp'~yhp'; call export(x,fname,names); */ endif; if plgr2; _ptek="c:\\archivos\\term_str\\svensson\\col\\colfig2.tkf"; graphprt("-c=4") ; title(" "\ "\LObserved and Estimated Bond Prices, Colombia, " $+ tdstring); if errbar; title(" "\ "\LObserved and Estimated Bond Prices, Colombia, " $+ tdstring $+ "\L95% confidence interval" ); endif; xlabel("Maturity"); ylabel(" "); _perrbar = 0; if errbarp; if nsx; if restr; b2 = bh[1 3 4 5 6]; yp = CONFNF(&PVNS22x,b2,cov,q,tol); else; yp = CONFNF(&PVNS2x,bh,cov,q,tol); endif; else; if restr; b2 = bh[1 3 4]; yp = CONFNF(&PVNS22,b2,cov,q,tol); else; b2 = bh[1:4]; yp = CONFNF(&PVNS2,b2,cov,q,tol); endif; endif; _perrbar = ttm+(tm'~tm'~tm')~yp[.,2 1 3]~ones(rows(tm'),1)*(6~0~0); endif; _plctrl = {-1,-1,-1,-1,1}; _pltype = {4,6,3,3,6,5}; _pstype = {2,8,4,5,4,6,7}; _psymsiz = {3.5,1,3.5,4,2,2,2}; _plegctl = {2,4,6.3,1}; _plegstr = "Observed bond price"\ "\000Estimated bond price"\ "\000Coupon rate + 100"; _pline = 0; scale(ttm+tm',pdat~ph~(100+c')); if xtics1 == 1; xtics(2000,2009,2,2); endif; if ytics1 == 1; @ytics(90,120,10,1);@ ytics(90,130,10,1); endif; _pmcolor={"","","","","","","","",15}; xy(ttm+tm',pdat~ph~(100+c')); endif; if output_file; /* Prices */ /* fname = "c:\\archivos\\term_str\\svensson\\col\\ts_p.xls"; let names = maturity oprice eprice; x = (ttm+tm')~pdat~ph; call export(x,fname,names); */ endif; if plgr3; title(" "\ "\LSpot and Forward Rates, Colombia, " $+ tdstring); if errbar; title(" "\ "\LSpot and Forward Rates, Colombia, " $+ tdstring $+ "\L95% confidence interval" ); endif; xlabel("Maturity"); ylabel("Annually compounded rate, %/yr"); _plctrl = {0,0,0,0,0,0}; if asy; _pline = 1~1~ttm+5~sinfh~ttm+40~sinfh~1~2~0; endif; _pcolor = {5,12}; _pltype = {3,6,6,1,5}; _pstype = {2,3,1,5,4,6,7}; _plegctl = {2,4,5.3,1}; _plegstr = "Estimated spot rate"\ "\000Estimated forward rate"; if errbar; if nsx; if restr; b2 = bh[1 3 4 5 6]; ys = 100*(exp(CONFNF(&spotr2x,b2,cov,q,tol)/100)-1); yf = 100*(exp(CONFNF(&forwr2x,b2,cov,q,tol)/100)-1); else; ys = 100*(exp(CONFNF(&spotr3x,bh,cov,q,tol)/100)-1); yf = 100*(exp(CONFNF(&forwr3x,bh,cov,q,tol)/100)-1); endif; else; if restr; b2 = bh[1 3 4]; ys = 100*(exp(CONFNF(&spotr2,b2,cov,q,tol)/100)-1); yf = 100*(exp(CONFNF(&forwr2,b2,cov,q,tol)/100)-1); else; b2 = bh[1:4]; ys = 100*(exp(CONFNF(&spotr3,b2,cov,q,tol)/100)-1); yf = 100*(exp(CONFNF(&forwr3,b2,cov,q,tol)/100)-1); endif; endif; _perrbar = ttm+(tm2~tm2~tm2)~ys[.,2 1 3]~ones(rows(tm2),1)*(6~0~0) | ttm+(tm2~tm2~tm2)~yf[.,2 1 3]~ones(rows(tm2),1)*(6~0~0) ; endif; scale(ttm+(tm1'~tm1'),(sh1'~fh1')); if xtics1 == 1; xtics(2000,2009,2,2); endif; if ytics1 == 1; ytics(10,23,1,1); endif; _pmcolor={"","","","","","","","",15}; xy(ttm+(tm1'~tm1'),(sh1'~fh1')); _pcolor = {5,3,12,7}; endif; if plgr4; if Discount_ci; @ ********************************************** @; @ Confidence intervals for the Discount Function @; @****** NS with Restriction ****@; proc SPOTR21(b2); local bb,s; bb = b2[1]|_s0-b2[1]|b2[2 3]; s = SPOTR(tm1',bb); retp(s); endp; proc DISCFN1(b); external proc SPOTR21; local s,d; s = SPOTR21(b); d = exp( -(s/100).*(tm1') ); retp( 100*d ); endp; @****** NS w/o Restriction ****@; proc SPOTR31(b); local s; s = SPOTR(tm1',b); retp(s); endp; proc DISCFN31(b); external proc SPOTR31; local s,d; s = SPOTR31(b); d = exp( -(s/100).*(tm1') ); retp( 100*d ); endp; @****** NSX with Restriction ****@; proc SPOTR2X1(b2); local bb,s; bb = b2[1]|_s0-b2[1]|b2[2 3 4 5]; s = SPOTRX(tm1',bb); retp(s); endp; proc DISCFNX1(b); external proc SPOTR2X1; local s,d; s = SPOTR2X1(b); d = exp( -(s/100).*(tm1') ); retp( 100*d ); endp; @****** NSX w/o Restriction ****@; proc SPOTR3X1(b); local s; s = SPOTRX(tm1',b); retp(s); endp; proc DISCFN3X1(b); external proc SPOTR3X1; local s,d; s = SPOTR3X1(b); d = exp( -(s/100).*(tm1') ); retp( 100*d ); endp; @******************************@; if nsx; if restr; b2 = bh[1 3 4 5 6]; dicfn_ci = confnf(&DISCFNX1,b2,cov,q,tol); else; dicfn_ci = confnf(&DISCFN3X1,b2,cov,q,tol); endif; else; if restr; b2 = bh[1 3 4]; dicfn_ci = confnf(&DISCFN1,b2,cov,q,tol); else; b2 = bh[1:4]; dicfn_ci = confnf(&DISCFN31,b2,cov,q,tol); endif; endif; @ ********************************************** @; title(" "\ "\LDiscount Function, Colombia, " $+ tdstring $+ "\L95% confidence interval"); xlabel("Maturity, yrs"); ylabel("%"); _perrbar = 0; _plctrl = {0,0,0}; _pcolor = {5,12,5}; _pltype = {3,6,3}; @_pstype = {9,8,8};@ _pline = 0; _plegctl = {2,4,1.5,1}; _plegstr = "Confidence Intervals"\ "\000Estimated discount function"\ "\000Confidence Intervals"; scale( (tm1'~ tm1'~ tm1') , dicfn_ci); /* xtics(0,17,1,1); */ ytics(40,110,10,1); _pmcolor={"","","","","","","","",15}; xy( (tm1'~ tm1'~ tm1') , dicfn_ci ); else; title(" "\ "\LDiscount Function, Colombia, " $+ tdstring); xlabel("Maturity, yrs"); ylabel("%"); _perrbar = 0; _plctrl = {0,-1,0,1}; _pltype = {6,6,3,6,1,5}; _pstype = {2,3,1,5,4,6,7}; _pline = 0; _plegctl = {2,4,2,1}; _plegstr = "Estimated discount function"; scale(tm1',ph1'); /* xtics(0,17,1,1); */ ytics(40,110,10,1); _pmcolor={"","","","","","","","",15}; xy(tm1',ph1'); endif; endif;