0001 function [mu, A, lambda_0, kappa, theta, sigma, rho] ...
0002 = heuristic_pcsv_param(y_t, dt, p)
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018 [T, n] = size(y_t);
0019
0020
0021 [A, Lambda, V] = svd(cov(diff(y_t))/dt);
0022 lambda_0 = diag(Lambda);
0023
0024 lambda_0 = lambda_0(1:p)';
0025 A = A(:,1:p);
0026
0027
0028 mu = (mean(diff(y_t))/dt + (0.5*(A.^2)*lambda_0')')';
0029
0030
0031 roll_period = 8;
0032 l = zeros(T-roll_period, p);
0033
0034 r_t = diff(y_t);
0035
0036
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
0043
0044
0045 end
0046
0047
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
0060 sigma(k) = sqrt(s/dt);
0061
0062 theta(k) = lambda_0(k);
0063 kappa(k) = -b(2);
0064
0065
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
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