0001 function [mu, Sigma_0, M, Q, rho, beta] = heuristic_wasc_param_2d(y_t, dt)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015 [T, n] = size(y_t);
0016 r = diff(y_t);
0017
0018 Sigma_0 = cov(r)/dt;
0019 mu = mean(r)'/dt + (0.5*diag(Sigma_0));
0020
0021
0022 roll_period = 25;
0023 S = zeros(n, n, T-roll_period);
0024 for k=1:(T-roll_period)
0025 S(:,:,k) = cov(diff(y_t(k:(k+roll_period-1),:)))/dt;
0026 end
0027
0028 dS = diff(S, 1, 3);
0029 K = size(dS,3);
0030 E = @(a) a'*mean(dS,3)*a;
0031 V = @(a) mean( (arrayfun(@(k) a'*dS(:,:,k)*a, 1:K) - E(a)).^2 );
0032 Cov = @(a,g) mean( ...
0033 (arrayfun(@(k) a'*dS(:,:,k)*a, 1:K) - E(a)) ...
0034 .* (arrayfun(@(k) g'*dS(:,:,k)*g, 1:K) - E(g)) );
0035 KS = size(S,3);
0036 ES = @(a) a'*mean(S,3)*a;
0037 VS = @(a) mean( (arrayfun(@(k) a'*S(:,:,k)*a, 1:KS) - ES(a)).^2 );
0038 K = @(a) 2*ES(a)^2 / VS(a);
0039 a1 = [1;0];
0040 a2 = [0;1];
0041 a12 = [1;1];
0042
0043 beta = K(a12);
0044 beta = max(beta, 9);
0045
0046 Q = zeros(n);
0047 Q(1,1) = sqrt( V(a1) / (4*Sigma_0(1,1)*dt) );
0048 Q(2,2) = sqrt( V(a2) / (4*Sigma_0(2,2)*dt) );
0049 Q(1,2) = -Cov(a1,a2) / (4*Sigma_0(1,2)*Q(1,1)*dt);
0050
0051 A = zeros(3);
0052 A(1,1) = 2*Sigma_0(1,1);
0053 A(1,2) = 2*Sigma_0(1,2);
0054 A(2,2) = 2*Sigma_0(1,2);
0055 A(2,3) = 2*Sigma_0(2,2);
0056 A(3,1) = 2*(Sigma_0(1,1)+Sigma_0(1,2));
0057 A(3,2) = 2*(Sigma_0(1,1)+2*Sigma_0(1,2)+Sigma_0(2,2));
0058 A(3,3) = 2*(Sigma_0(2,2)+Sigma_0(1,2));
0059 b = zeros(3,1);
0060 b(1) = E(a1)/dt - beta*a1'*(Q'*Q)*a1;
0061 b(2) = E(a2)/dt - beta*a2'*(Q'*Q)*a2;
0062 b(3) = E(a12)/dt - beta*a12'*(Q'*Q)*a12;
0063 m = A \ b;
0064 M = [m(1), m(2); m(2), m(3)];
0065
0066 period_length = 50;
0067 number_of_periods = floor((T-1)/period_length)-1;
0068 dZ = zeros(number_of_periods*n, 1);
0069 dW = zeros(number_of_periods*n, n);
0070 S_last = Sigma_0;
0071 for k=1:number_of_periods
0072 idx = (k*period_length):((k+1)*period_length);
0073 period_returns = r(idx, :);
0074 S = cov(period_returns);
0075 dS = (S - S_last) / period_length;
0076 dY = mean(period_returns);
0077 sqrt_S = chol(S);
0078 dZ((n*(k-1)+1):(n*k)) = sqrt_S \ (dY' - (mu - diag(S))*dt);
0079 dW((n*(k-1)+1):(n*k), :) = (2*sqrt_S*Q) \ (dS-(beta*(Q'*Q)+2*M*S)*dt);
0080 end
0081 rho = (dW \ dZ);
0082 end