Home > code > core > estimation > heuristic_pcsv_param.m

heuristic_pcsv_param

PURPOSE ^

HEURISTIC_PCSV_PARAM Obtains heuristic parameter estimates for a PCSV model.

SYNOPSIS ^

function [mu, A, lambda_0, kappa, theta, sigma, rho]= heuristic_pcsv_param(y_t, dt, p)

DESCRIPTION ^

HEURISTIC_PCSV_PARAM Obtains heuristic parameter estimates for a PCSV model.

  [mu, A, lambda_0, kappa, theta, sigma, rho]
    = heuristic_pcsv_param(y_t, dt) calculates parameter estimates for a PCSV
    model with the time series y_t using a heuristic approach. The results
    are supposed to be used as starting points for an iterative optimization
    algorithm.

    INPUT y_t: A Txn matrix representing the log prices
           dt: The time difference in the observed series y_t
            p: Number of eigenvectors used for this model

 created by Benedikt Rudolph
 DATE: 16-Aug-2012

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [mu, A, lambda_0, kappa, theta, sigma, rho] ...
0002           = heuristic_pcsv_param(y_t, dt, p)
0003 %HEURISTIC_PCSV_PARAM Obtains heuristic parameter estimates for a PCSV model.
0004 %
0005 %  [mu, A, lambda_0, kappa, theta, sigma, rho]
0006 %    = heuristic_pcsv_param(y_t, dt) calculates parameter estimates for a PCSV
0007 %    model with the time series y_t using a heuristic approach. The results
0008 %    are supposed to be used as starting points for an iterative optimization
0009 %    algorithm.
0010 %
0011 %    INPUT y_t: A Txn matrix representing the log prices
0012 %           dt: The time difference in the observed series y_t
0013 %            p: Number of eigenvectors used for this model
0014 %
0015 % created by Benedikt Rudolph
0016 % DATE: 16-Aug-2012
0017 
0018   [T, n] = size(y_t);
0019   
0020   % obtain A and lambda_0
0021   [A, Lambda, V] = svd(cov(diff(y_t))/dt);
0022   lambda_0 = diag(Lambda);
0023   % truncate to desired number of eigenvectors/eigenvalues
0024   lambda_0 = lambda_0(1:p)';
0025   A = A(:,1:p);
0026   
0027   % obtain mu
0028   mu = (mean(diff(y_t))/dt + (0.5*(A.^2)*lambda_0')')';
0029   
0030   % obtain a time series l of eigenvalues using rolling estimates
0031   roll_period = 8;
0032   l = zeros(T-roll_period, p);
0033   %_dY = zeros(T-roll_period, p);
0034   r_t = diff(y_t);
0035   %dB = zeros(T-roll_period-1, p);
0036   %dW = zeros(T-roll_period-1, p);
0037   for k=1:(T-roll_period)
0038     S = cov(diff(y_t(k:(k+roll_period-1),:)))/dt;
0039     [A_k, Lambda_k, V_k] = svd(S);
0040     lambda_k = diag(Lambda_k);
0041     l(k,:) = lambda_k(1:p);
0042     %_dY(k,:) = mean(diff(y_t(k:(k+roll_period-1),:)));
0043     %dW(k,:) = (A*diag(sqrt(l(k,:))))*sqrt(dt) ...
0044     %         \ (r_t(k,:)' - (mu - 0.5*A.^2*l(k,:)')*dt);
0045   end
0046   
0047   % Obtain kappa,theta and sigma using OLS on the eigenvalue estimates l
0048   kappa = zeros(1,p);
0049   theta = zeros(1,p);
0050   sigma = zeros(1,p);
0051   feller_condition = @(kappa, theta, sigma) 2*kappa.*theta > sigma.^2;
0052   for k=1:p
0053     lambda = l(:,k);
0054     Y = diff(lambda)./sqrt(lambda(1:(end-1)));
0055     X = [dt./sqrt(lambda(1:(end-1))), dt*sqrt(lambda(1:(end-1)))];
0056     b = X \ Y;
0057     r = Y - X*b;
0058     s = (r'*r) / (size(X,1)-2);
0059     %dB(:,k) = r/sqrt(s);
0060     sigma(k) = sqrt(s/dt);
0061     %theta(k) = -b(1) / b(2);
0062     theta(k) = lambda_0(k);
0063     kappa(k) = -b(2);
0064     % if kappa, theta and sigma do not satisfy the Feller condition
0065     % shift the parameters equally so that they do satisfy it
0066     if ~feller_condition(kappa(k), theta(k), sigma(k))
0067       alpha = - ( sqrt(2*kappa(k).*theta(k)) - sigma(k) ) ./ ...
0068                 ( sqrt(2*kappa(k).*theta(k)) + sigma(k) ) * 1.01;
0069       kappa(k) = (1+alpha) * kappa(k);
0070       theta(k) = (1+alpha) * theta(k);
0071       sigma(k) = (1-alpha) * sigma(k);
0072     end
0073   end
0074   
0075    % Obtain rho
0076   period_length = 12;
0077   number_of_periods = floor((T-1)/period_length)-1;
0078   dB = zeros(number_of_periods, p);
0079   dW = zeros(number_of_periods, p);
0080   lambda_last = lambda_0;
0081   for k=1:number_of_periods
0082     idx = (k*period_length):((k+1)*period_length);
0083     period_returns = r_t(idx, :);
0084     S = cov(period_returns);
0085     [A_k, Lambda_k, V_k] = svd(S);
0086     lambda_k = diag(Lambda_k);
0087     lambda_k = lambda_k(1:p);
0088     dlambda = (lambda_k'-lambda_last) / period_length;
0089     dY = mean(period_returns);
0090     dB(k,:) = ( (dlambda - kappa.*(theta-lambda_k')*dt) ./ ...
0091                                   (sigma .* sqrt(lambda_k')) ) / sqrt(dt);
0092     dW(k, :) = ( (A*diag(sqrt(lambda_k)))*sqrt(dt) ...
0093               \ (dY' - (mu - 0.5*A.^2*lambda_k)*dt) )' / sqrt(dt);
0094   end
0095   rho = diag(corr(dW,dB))';
0096 end

Generated on Mon 29-Apr-2013 19:29:13 by m2html © 2005