如何使用MATLAB进行分位数回归？

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

各位童鞋：

        本人正在研究分位数回归在金融领域的应用，看了好多资料上的分位数回归都是采用R、SAS等其他统计软件进行分位数回归的参数估计，请问使用MATLAB如何进行分位数回归？MATLAB中有没有现成的分位数估计程序？对分位数回归感兴趣的同志可以加QQ393419199一起讨论，谢谢！

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏11 回帖

关键词：MATLAB 分位数回归 matla atlab 如何使用 如何

见过这个方法，可是自己不会用。有什么书籍介绍这个方法吗？

楼主自己看看把

liwenxue_137 发表于 2013-10-15 10:33
楼主自己看看把

写错了，这个只是分位数算法，稍微修改就可以了

liwenxue_137 发表于 2013-10-15 10:33
楼主自己看看把

谢谢你的M文件，不过这个M文件是求向量的分位数的，不是分位数回归，MATLAB应该没有现成的分位数回归程序。

没有查过，我用的是R软件做的，有个quantreg分位数包，应有尽有

我前一个月也用matlab，现在发现matlab做统计，很多包包都没有，就该R了。学习R只需要几天就可以了

liwenxue_137 发表于 2013-10-30 08:22
我前一个月也用matlab，现在发现matlab做统计，很多包包都没有，就该R了。学习R只需要几天就可以了

谢谢，stata也可以做分位数回归，我不是很熟悉R

function b = rq_fnm(X, y, p)
% Construct the dual problem of quantile regression
% Solve it with lp_fnm
%
%
[m n] = size(X);
u = ones(m, 1);
a = (1 - p) .* u;
b = -lp_fnm(X', -y', X' * a, u, a)';

function y = lp_fnm(A, c, b, u, x)
% Solve a linear program by the interior point method:
% min(c * u), s.t. A * x = b and 0 < x < u
% An initial feasible solution has to be provided as x
%
% Function lp_fnm of Daniel Morillo & Roger Koenker
% Translated from Ox to Matlab by Paul Eilers 1999
% Modified by Roger Koenker 2000--
% More changes by Paul Eilers 2004


% Set some constants
beta = 0.9995;
small = 1e-5;
max_it = 50;
[m n] = size(A);

% Generate inital feasible point
s = u - x;
y = (A' \  c')';
r = c - y * A;
r = r + 0.001 * (r == 0);    % PE 2004
z = r .* (r > 0);
w = z - r;
gap = c * x - y * b + w * u;

% Start iterations
it = 0;
while (gap) > small & it < max_it
    it = it + 1;
   
    %   Compute affine step
    q = 1 ./ (z' ./ x + w' ./ s);
    r = z - w;
    Q = spdiags(sqrt(q), 0, n, n);
    AQ = A * Q;          % PE 2004
    rhs = Q * r';        % "
    dy = (AQ' \ rhs)';   % "
    dx = q .* (dy * A - r)';
    ds = -dx;
    dz = -z .* (1 + dx ./ x)';
    dw = -w .* (1 + ds ./ s)';
   
    % Compute maximum allowable step lengths
    fx = bound(x, dx);
    fs = bound(s, ds);
    fw = bound(w, dw);
    fz = bound(z, dz);
    fp = min(fx, fs);
    fd = min(fw, fz);
    fp = min(min(beta * fp), 1);
    fd = min(min(beta * fd), 1);
   
    % If full step is feasible, take it. Otherwise modify it
    if min(fp, fd) < 1
        
        % Update mu
        mu = z * x + w * s;
        g = (z + fd * dz) * (x + fp * dx) + (w + fd * dw) * (s + fp * ds);
        mu = mu * (g / mu) ^3 / ( 2* n);
        
        % Compute modified step
        dxdz = dx .* dz';
        dsdw = ds .* dw';
        xinv = 1 ./ x;
        sinv = 1 ./ s;
        xi = mu * (xinv - sinv);
        rhs = rhs + Q * (dxdz - dsdw - xi);
        dy = (AQ' \ rhs)';
        dx = q .* (A' * dy' + xi - r' -dxdz + dsdw);
        ds = -dx;
        dz = mu * xinv' - z - xinv' .* z .* dx' - dxdz';
        dw = mu * sinv' - w - sinv' .* w .* ds' - dsdw';
        
        % Compute maximum allowable step lengths
        fx = bound(x, dx);
        fs = bound(s, ds);
        fw = bound(w, dw);
        fz = bound(z, dz);
        fp = min(fx, fs);
        fd = min(fw, fz);
        fp = min(min(beta * fp), 1);
        fd = min(min(beta * fd), 1);
        
    end
   
    % Take the step
    x = x + fp * dx;
    s = s + fp * ds;
    y = y + fd * dy;
    w = w + fd * dw;
    z = z + fd * dz;
    gap = c * x - y * b + w * u;
    %disp(gap);
end

function b = bound(x, dx)
% Fill vector with allowed step lengths
% Support function for lp_fnm
b = 1e20 + 0 * x;
f = find(dx < 0);
b(f) = -x(f) ./ dx(f);

这个就是matlab做分位数回归。最近在看这个，R，matlab都看了个够，2个月了

function [p,stats]=quantreg(x,y,tau,order,Nboot)
% Quantile Regression
%
% USAGE: [p,stats]=quantreg(x,y,tau[,order,nboot]);
%
% INPUTS:
%   x,y: data that is fitted. (x and y should be columns)
%        Note: that if x is a matrix with several columns then multiple
%        linear regression is used and the "order" argument is not used.
%   tau: quantile used in regression.
%   order: polynomial order. (default=1)
%          (negative orders are interpreted as zero intercept)
%   nboot: number of bootstrap surrogates used in statistical inference.(default=200)
%
% stats is a structure with the following fields:
%      .pse:    standard error on p. (not independent)
%      .pboot:  the bootstrapped polynomial coefficients.
%      .yfitci: 95% confidence interval on "polyval(p,x)" or "x*p"
%
% [If no output arguments are specified then the code will attempt to make a default test figure
% for convenience, which may not be appropriate for your data (especially if x is not sorted).]
%
% Note: uses bootstrap on residuals for statistical inference. (see help bootstrp)
% check also: http://www.econ.uiuc.edu/~roger/research/intro/rq.pdf
%
% EXAMPLE:
% x=(1:1000)';
% y=randn(size(x)).*(1+x/300)+(x/300).^2;
% [p,stats]=quantreg(x,y,.9,2);
% plot(x,y,x,polyval(p,x),x,stats.yfitci,'k:')
% legend('data','2nd order 90th percentile fit','95% confidence interval','location','best')
%
% For references on the method check e.g. and refs therein:
% http://www.econ.uiuc.edu/~roger/research/rq/QRJEP.pdf
%
%  Copyright (C) 2008, Aslak Grinsted

%   This software may be used, copied, or redistributed as long as it is not
%   sold and this copyright notice is reproduced on each copy made.  This
%   routine is provided as is without any express or implied warranties
%   whatsoever.


if nargin<3
    error('Not enough input arguments.');
end
if nargin<4, order=[]; end
if nargin<5, Nboot=200; end

if (tau<=0)||(tau>=1),
    error('the percentile (tau) must be between 0 and 1.')
end

if size(x,1)~=size(y,1)
    error('length of x and y must be the same.');
end

if numel(y)~=size(y,1)
    error('y must be a column vector.')
end

if size(x,2)==1
    if isempty(order)
        order=1;
    end
    %Construct Vandermonde matrix.
    if order>0
        x(:,order+1)=1;
    else
        order=abs(order);
    end
    x(:,order)=x(:,1); %flipped LR instead of
    for ii=order-1:-1:1
        x(:,ii)=x(:,order).*x(:,ii+1);
    end
elseif isempty(order)
    order=1; %used for default plotting
else
    error('Can not use multi-column x and specify an order argument.');
end


pmean=x\y; %Start with an OLS regression

rho=@(r)sum(abs(r-(r<0).*r/tau));

p=fminsearch(@(p)rho(y-x*p),pmean);

if nargout==0
    [xx,six]=sortrows(x(:,order));
    plot(xx,y(six),'.',x(six,order),x(six,:)*p,'r.-')
    legend('data',sprintf('quantreg-fit ptile=%.0f%%',tau*100),'location','best')
    clear p
    return
end

if nargout>1
    %calculate confidence intervals using bootstrap on residuals
   
    yfit=x*p;
    resid=y-yfit;
   
    stats.pboot=bootstrp(Nboot,@(bootr)fminsearch(@(p)rho(yfit+bootr-x*p),p)', resid);
    stats.pse=std(stats.pboot);
   
    qq=zeros(size(x,1),Nboot);
    for ii=1:Nboot
        qq(:,ii)=x*stats.pboot(ii,:)';
    end
    stats.yfitci=prctile(qq',[2.5 97.5])';
   
end


%
%
%
% function ps=invtranspoly(p,kx)
%
% ps=p;
% order=length(p);
% fact=cumprod(1:order);
% kx=cumprod(repmat(kx,1,order));
% for ii=2:order
%     for jj=1:ii-1
%         N=order-jj+1;
%         k=ii-jj;
%         k=min(k,N-k);
%         ps(ii)=ps(ii)+p(jj)*kk(k)*fact(N)/(fact(k)*fact(N-k));
%     end
% end