基于SIC准则的重庆市私人汽车的结构变化模型

基于SIC准则的重庆市私人汽车的结构变化模型

重庆大学 郭雪梅、明月、马艳娜

摘要：重庆作为共和国最年轻的、西部地区唯一的直辖市，具有承东启西的地区优势，研究重庆经济结构变化有着重大的意义。本文以重庆统计年鉴私人汽车拥有量15年的数据为研究材料，通过回归分析及主成分回归模型得出了关于重庆市私人汽车拥有量的回归方程，该方程表明：公路里程、城市居民家庭人均可支配收、重庆市生产总值、重庆市总人口对重庆市私人汽车拥有量的贡献系数均为正，其中公路里程的方差贡献率最大。再由Schwarz信息准则方法准确找出了重庆市汽车产业结构在2002年发生了变化，并从各方面分析了造成该点发生转变的原因。最后对重庆市汽车产业的发展提出建议。

关键词：汽车拥有量 回归模型 主成分分析 转变点 Schwarz信息准则

一、引言

（一）重庆市汽车产业的现状及发展

汽车产业的发展水平在一定程度上是衡量一个国家经济实力、工业水平和科技创新能力高低的重要标志之一，其对国民经济和社会发展的巨大推动作用使得许多工业发达地区都将其作为支柱产业并加以重点发展和扶植。汽车产业作为重庆市的第一支柱产业，经过40多年的发展，目前，重庆已拥有汽车企业24家，行业资产超过600亿元。自“十一五”以来，重庆汽车产业保持了良好的发展势头，产量位居全国前列。重庆汽车产业对推动重庆区域经济的发展将产生巨大作用。

2010年重庆市实现地区生产总值7894.24亿元，比上年增长17.1%。其中，第一产业实现增加值685.39亿元，增长6.1%。私人汽车保有量77.03万辆，增长33.7%，私人轿车保有量39.93万辆，增长36.3%。全年新注册汽车27.55万辆，增长33.8%；其中，新注册轿车13.86万辆，增长33.0%。2010年，重庆全年生产汽车220万辆，位列全国汽车产量前三。车市销售突破600亿元，增幅继续领跑。预计未来5-10年，重庆汽车销量增长将保持20%以上。

（二）基于重庆市汽车产业目前的研究情况

目前国内学者更多地关注重庆汽车产业经营过程中的市场价格变化、市场预测、销售情况、汽车行业的发展策略、汽车产量增长性分析及相关政策法规对行业的影响等问题。比如：毕建平、王新莉、邵凤翔、马建辉对中国民用汽车拥有量进行时间序列动态模型的研究；徐小钦、杨红艳关于重庆市汽车产业集群发展的对策研究等。

（三）本文模型的构架

本文通过建立Change point模型架构,分析转变点时期的重大事件对汽车产业发展以及对重庆市经济结构变化的影响。由于各变量之间往往存在着一定程度的相关性，而人们希望能从这些分析指标中尽可能多地提取有效信息。本文用因子分析、主成分分析的方法通过R软件运行计算研究影响民用私有汽车拥有量的因素以及用Matlab软件进行转变点的寻找。在研究过程中讨论各变量对重庆汽车拥有量的相关性。本文收集了1995年到2009年15年的数据，包括:汽车拥有量、居民可支配收入、GDP、公路里程，通过R软件得出回归模型，利用转变点分析法找出影响汽车拥有量的主要因素。最后通过我们的模型分析重庆市私有汽车消费现状的提出发展汽车行业的政策建议。

二、资料来源及描述

（一）资料来源背景

重庆经济发展迅速，汽车行业尤为突出，私人汽车拥有量逐年增加，但汽车不属于生活必需品的消费，因此购买汽车的费用就与家庭可支配收入直接关联，其次，汽车产业的发展更多的是来自国家、政府宏观经济政策的扶持，因此GDP也是一个有效衡量汽车产业发展的重要指标。与此同时，人们出行越来越多地依赖汽车，对公路里程、公路质量的要求逐步加高。现在重庆各区县都积极修建高速公路、城市道路、乡村道路、桥梁及隧道。选择公路里程指标来说明私人汽车拥有量的增长与道路有关,公路长度的增加也在影响着汽车的销量。

（二）数据及数据的基本特征

从表1可以看出，在1995年-2009年间，重庆市总人口、城市居民家庭人均可支配收入与重庆市生产总值增长平稳，没有出现较大波动。但公路里程与民用汽车拥有量在某年份跳跃幅度比较大。

在收集数据时，能找到的关于私人民用汽车拥有量的最早数据是1995年，所以样本数据量较小。除了上述因素外，还有一些重要的影响因素，比如说品牌及国家政策对汽车拥有量的影响。由于样本数据量有限，不能选择过多的解释变量，以避免造成模型自由度为负的错误。

本文数据均来自《重庆市统计年鉴》。

表1 样本原始数据表

其中表1为原始数据，表2中、、、、是对表1中、、、、等数据取以10为底的对数。

表2 样本处理数据表

三、模型的建立及求解

（一）回归分析

回归分析(regression analysis)是研究一个变量关于另一个（些）变量的具体依赖关系的计算方法和理论， 从一组样本数据出发，确定变量之间的数学关系式对这些关系式的可信程度进行各种统计检验，并从影响某一特定变量的诸多变量中找出变量的显著性，利用所求的关系式，根据一个或几个变量的取值来预测或控制另一个特定变量的取值，并给出这种预测或控制的精确程度。

根据表2的数据，建立回归模型，分析公路线路里程（公里）、城市居民家庭人均可支配收入（元）、重庆市生产总值（亿元）、重庆市总人口（万人）对重庆市私人车辆拥有量（辆）的影响。建立回归模型如下：

假定

其中表示随、、、的变化而线性变化的部分；为随机误差，它是其他一切不确定因素的总和，其值不可观测，通常假定。称、、、为回归自变量，为回归因变量，为回归常数，、、、为回归系数。

使用R软件对表二的数据进行回归分析并分别画出、、、与的关系图如下：

图1 各变量与私人车拥有量关系图

先做表二数据的回归分析，得到下面的输出：

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 188.65723 182.24760 1.035 0.325

x1 -0.13805 0.58320 -0.237 0.818

x2 -0.03798 6.93596 -0.005 0.996

x3 5.87879 5.59459 1.051 0.318

x4 -57.95460 55.03946 -1.053 0.317

从输出的结果可以看出，按4个变量得到的回归方程为：

（3.1.1）

分析式（3.1.1），发现它并不合理，回到问题本身，是私人车辆拥有量，是公路线路里程，对应系数的符号为负，即公路线路里程越多，私人车辆拥有量越少，这与实际情况不符；是城市居民家庭人均可支配收入，对应系数的符号为负，也就是说城市居民家庭人均可支配收入越多，私人车辆拥有量越少，这与实际情况相悖。问其原因，4个变量存在着多重共线性，进而对数据进行主成分分析，找出其变量间的共线性关系

（二）主成分分析

主成分分析（Principal Component Analysis，PCA）是将多个变量通过线性变换以选出较少个数重要变量的一种多元统计分析方法，这些主成分能够反映原始变量的绝大部分信息，它们通常表示为原始变量的某种线性组合。其基本思想是设法将原来众多具有一定相关性的指标重新组合成一组新的相互独立或不相关的综合指标来代替原来的指标。通常数学上的处理就是将原来的指标作线性组合，作为新的综合指标。具体步骤为：1. 指标数据标准化（R软件自动执行）；2. 指标之间的相关性判定；3. 确定主成分个数m；4.主成分表达式；5. 主成分命名。

对表二用R软件进行主成分分析结果如下：

Importance of components:

Comp.1 Comp.2 Comp.3 Comp.4

Standard deviation 1.968206 0.34596826 0.073701602 0.0322208910

Proportion of Variance 0.968459 0.02992351 0.001357982 0.0002595465

前三个主成分已达到99.97%的贡献率，而且=，所以变量存在着多重共线性 。

为克服多重共线性的影响，对变量做主成分回归:

(A)：预测样本主成分并做主成分分析，在R软件中运算结果如下

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 5.27851 0.05892 89.58 < 2e-16 ***

Z 0.32973 0.02994 11.01 5.8e-08 ***

Signif. codes: 0 ‘***’0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’0.1 ‘ ’1

Residual standard error: 0.2282 on 13 degrees of freedom

Multiple R-squared: 0.9032, Adjusted R-squared: 0.8958

F-statistic: 121.3 on 1 and 13 DF, p-value: 5.802e-08

回归系数和回归方程均通过检验，而且效果显著，即得到回归方程为

（3.2.1）

上述方程得到的是相应变量与主成分的关系，但还是希望得到相应变量与原变量之间的关系。

(B)作变换，得到原坐标下的关系表达式在R软件中运算结果如下

(Intercept) 　　 x1 　 x2 　 x3 x4

-54.6081213 　　 0.5962799 　　 0.9116767 　　 0.7234347 　 14.6212501

即得到回归方程为

（3.2.2）

此时相应的因变量系数均为正数，比原回归方程（3.1.1）更合理。

（三）转变点分析

由上述回归方程（3.2.2）可以得到各自变量与因变量呈线性关系，然而从图一中可以看出各变量之间并非有较好的线性关系，而是在某段时间跳跃性很大，因此我们需要研究在哪些点发生了实质性的变化。其中Schwarz信息准则方法（SIC）是一种有效找出各转变点的方法，它用于检测转变点的存在性时，无需导出其复杂的分布函数，所以利用信息准则来估计变点的个数和位置是较为简单的。

Schwarz信息准则定义为：， 这里是极大似然函数，是是参数的极大似然估计，p是模型中的自由参数的个数，n是样本大小。按此定义，在假设H0成立下的SIC记为SIC(n)，而在H1成立下的SIC记为SIC(k)，选择接受H0或H1将取决于最小信息准则的原理。

为了描述在其它条件不变的情况下，私人民用车辆拥有量对于某个时间发生扰动的敏感程度，我们引入以下模型：

（3.3.1）

其中为白噪声，为对的相依水平。模型（3.3.1）的转变点问题可以表述为如下的假设检验：

其中是带估的转变点位置， =，我们需要判断从1到m的期间内私人汽车拥有量结构是否发生了变化，并在哪个时段发生了变化。我们根据SCI来估计转变点的位置，只需要估算出值与及即可。

记下的SIC为SIC（m），则：

其中

,

记下的SIC为SIC(k)，有

其中　

采用SIC准则进行判决：

（1）若，则不拒绝，即认为没有证据表明到时间m有转变发生；

（2）若，则拒绝，即认为在时间发生转变，其中满足．则是真实转变点的一致估计．

对表2有如下假设：

原假设

备择假设

通过MATLAB软件计算有，和各个得出的结果如表3：

表3 SIC计算结果表

又， ，， 即2002年是转变点。

我们分别对转变前和转变后的数据做回归仍然有上述回归分析、主成分分析和主成分回归，运行结果如下：

转变前回归分析结果如下：

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) -41.0488 64.9499 -0.632 0.572

x1 -2.0228 1.1787 -1.716 0.185

x2 2.3097 2.1978 1.051 0.370

x3 0.8863 1.6508 0.537 0.629

x4 12.3888 19.8940 0.623 0.578

主成分分析结果：

Importance of components:

Comp.1 Comp.2 Comp.3 Comp.4

Standard deviation 1.9742576 0.28382974 0.124364807 0.079253405

Proportion of Variance 0.9744232 0.02013983 0.003866651 0.001570276

Cumulative Proportion 0.9744232 0.99456307 0.998429724 1.000000000

公路里程仍然有第一主成分，居民人均收入，GDP和人口之间有较大的相关性。

主成分回归结果：

(Intercept) x1 x2 x3 x4

-41.4967796 1.2204327 0.8179838 0.5931995 10.2585662

所以转变前的回归方程为：

转变前个变量与汽车拥有量的关系图如下：

图2 转变前各变量与民用车拥有量关系图

转变后回归分析结果如下：

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) -59.0151 115.1263 -0.513 0.659

x1 0.1870 0.1032 1.812 0.212

x2 -2.3230 4.7045 -0.494 0.670

x3 1.7886 2.3208 0.771 0.521

x4 19.1398 36.0934 0.530 0.649

这样的效果并不理想，究其原因是变量之间有较大的相关性。

主成分分析：

Importance of components:

Comp.1 Comp.2 Comp.3 Comp.4

Standard deviation 1.929886 0.51715404 0.087519626 0.020791440

Proportion of Variance 0.931115 0.06686208 0.001914921 0.000108071

Cumulative Proportion 0.931115 0.99797701 0.999891929 1.000000000

仍然有公路里程是第一主成分，居民人均收入，GDP和人口之间有较大的相关性。

做主成分回归结果如下：

(Intercept) x1 x2 x3 x4

-17.9675152 0.1599558 0.3701295 0.3025412 5.8677580

所以转变后的回归方程为：

转变后个变量与汽车拥有量的关系图如下：

图3 转变后各变量与私人车拥有量关系图

可以对比得出结论图2、图3比图1的线性性更优，即转变点8前后经济结构发生了明显的变化。又由表3和图1可以看出2001年与2002年之间存在很大跳跃性，因此着重分析2001年及2002年的相关事件。

注：上述程序依次见附录1

四、转变点分析

（一）公路里程的影响

从上述回归模型：

（3.2.2）中可以看出，公路线路里程对私人汽车拥有量具有正的贡献系数，此系数为0.5962799，表明公路线路里程的不断增加对私有汽车拥有量的增加有比较大的推动作用。这也是符合实际情况的，公路增加了，人们可以利用较好的交通环境得到便捷，拥有汽车也会逐渐成为人们心里的想法，从而去购买它。我们做出了公路里程的增长率和私人车辆拥有量增长率如下表：

表4 增长率关系表

根据表4，做出下图：

图4 公路增长率与汽车增长率对比图

随着年数的增加，公路里程数不断增加，在2002年后，增长的比率都比较大，重庆市自2001年和2002 年起，加快了城市道路工程的修建。城市公路的修建对于城市经济发展以及民生工程建设都有很大的影响，单从上述主成分分析得到的数据就可以看出，公路里程对于私人汽车拥有量的增长有很大的贡献。其贡献率达到96.8%。因此公路里程可以作为汽车拥有量影响因素的主成分。在2001年后私人车辆拥有量增长幅度加大。汽车拥有量的增加也推动了重庆市的公路建设。自2001年颁布了一系列城市公路建设措施之后，汽车拥有量增加幅度稳定上升。

重庆以山城闻名，城市交通一直是发展的瓶颈。从2001年起，重庆市委、市政府颁布了一系列城市交通畅通工程实施政策如下：

（1）．《重庆市人民政府办公厅关于转发重庆市2001年城市道路交通畅通工程实施意见的通知》(渝办发[2001]23)

市委、市政府把实施城市道路交通畅通工程（以下简称“畅通工程”）列为2001年为民办实事八大“民心工程”和西部大开发10件大事之一。

（2）．《重庆市人民政府办公厅转发市公安局等部门关于重庆市2002年城市道路畅通工程实施意见补充方案的通知》

为认真贯彻落实重庆市第二次党代会精神，加快以“八小时重庆”、“半小时主城”交通为重点的基础设施建设，进一步巩固和深化实施畅通工程，全面推进富民兴渝，加快建设长江上游经济中心提供更加安全、畅通、有序、文明的道路交通秩序和环境，根据公安部、建设部《2002年全国城市道路交通管理畅通工程总体方案》要求，特制定重庆市2002年城市道路畅通工程实施意见的补充方案。

（3）．《重庆市人民政府关于认真做好2002年市重点建设项目和重大前期工作项目有关工作的通知》

交通工程共12项，总投资6052090万元，其中包括各区县路网、道桥及隧道建设工程。

（二）人均可支配收入的影响

随着经济的高速发展，居民收入水平提高，购买力增强，私有汽车已逐渐由高档奢侈品成为相对应的中高档奢侈品。人民可支配收入和生产总值的增长已成为私有汽车销量增加不可忽略的原因。从上述的模型 （3.2.2）可以看出人均可支配收入和生产总值对汽车拥有量的贡献分别是0.9116767和0.7234347。这是非常符合实际情况的，只有在人们有足够多的物质财富时，才会向往更好的生活方式。

下表给出了1995年到2009年人均收入增长率与私人车辆拥有量增长率的情况。

表5 各因素增长率关系表

根据表5，做出下图：

图5 重庆生产总值、人均可支配收入与汽车增长率

从图中可以看出，重庆市人均可支配收入每年基本呈现出均衡的增长速率，自中国2001年加入WTO，与国际接轨，加大了与国外的贸易交流，各地区的GDP也不断的增长，重庆从2001年开始，GDP以每年10%以上不断增长，人均可支配收入增长率与GDP的增长率基本相平。人民拥有足够的物质财富之后，自然会追求精神上的满足。

目前汽车产品在我国消费市场中属于高档耐用品，中低层家庭收入几乎都用于购买生活必须品以维持生计，无力购买像汽车这样的高档品，消费受到高度压抑。由于购买力不足，消费结构升级缓慢，收入差距两极分化现象日趋严重，这些对我国汽车消费市场产生负面影响。因此，改变中低家庭的消费结构，降低收入差距对提高私有汽车购买量产生间接影响，以刺激汽车产业的发展。

五、重庆市汽车消费市场现状及发展的制约因素

从上述的转变点分析可以看到，公路里程、人均收入、GDP对私人汽车的拥有量有比较大的推动作用，私人汽车拥有量在2001年、2002有较大的增加。反之，私人汽车拥有量也推动了城市建设，体现了一个城市的基础设施建设体系的完备状况；也部分体现了人民生活的富裕程度和经济发展状况。汽车产业带动了经济的增长，一个地区如果大力带动本土汽车产业发展，GDP、人均收入自然会增加，城市建设也会更加完备，人民也会越来越富裕。

虽然人们拥有的私有汽车量越来越多，但长安汽车、奇瑞、比亚迪、中华、吉利等我国自主汽车品牌的市场占有份额仅为25%。造成此种结果的原因首先是人们的消费观念倾向于国外汽车产品。其次是国内汽车品牌自主创新能力差，其产品主要服务于中低端市场，虽然国产汽车价格低廉、使用成本低，但其动力性能差、车身工艺粗糙、内饰材质普遍不够精致、安全性能不稳定；现如今由于收入水平的提高，思想观念的不同，使得人们更倾向于品质上的享受。再次，自中国加入WTO后，国内自主汽车品牌遭受国外汽车品牌的冲击使得其占有份额下降。最后，国外汽车企业历史悠久、工艺精湛、发展成熟、售后服务完善，这也是人们倾向于购买国外汽车品牌的原因之一。这种情况同样也适用于重庆市汽车消费市场。

目前重庆汽车较为成熟的车型多属于20世纪80年代引进的，如重型车引进法国贝利埃技术和奥地利斯泰尔技术，轻型车引进日本铃木技术，微型车同样引进日本铃木技术，像奥拓、羚羊等车型，在国内已生产十几年，与国内外同类车型的竞争中，技术上已明显落后。国外企业为保护其在中国推出的车型实现最大的利润率，通常会推迟引入新工艺或推广新产品。因此要想在技术上处于领先，根本上还要靠本地企业的自主创新。但目前，重庆汽车企业多数存在对国外技术的依赖，很多企业仍以短期获利以最大的全散件组装方式为主。

六、振兴重庆汽车产业的建议 虽然重庆市汽车产业发展日渐成熟，呈现规模化效应，但重庆将70%的投资都用于从国外进口发动机生产线、购买设备。因此，振兴重庆汽车产业我们提出一下建议：（一）加强自主创新、知识产权保护，高度重视专利工作 重庆汽车产业应在有市场针对性的前提下，加强与高校的科研合作，引进国内外的创新型人才，增强技术的专业性和可行性，提高自主创新能力，并且注重知识产权的保护。在不断吸收国外先进技术及管理销售方法时，要集中产业优势调高企业核心竞争力，加强自主研发创新的意识，并高度重视专利保护工作。

（二）准确定位产品细分市场，通过产品差异提高市场占有率

分清市场需求，有针对性的对于不同层次，不同年龄段的人们对不同类型汽车需求的供应。在提高产品质量的同时提高售后服务需求，在拓展新客户的同时，留住老客户。通过非对产品的细分，提高市场占有率，创造经济收入。

（三）通过企业并购、联合，形成规模经济

重庆作为我国六大老工业基地之一，许多国有制造业企业的机械工程配套能力强，具有完备的制造业基础，为重庆汽车产业的发展创造了良好的条件。因此，政府应给予汽车零部件企业集团政策性扶植，同时，鼓励国内外大型汽车企业对重庆市汽车零部件企业进行兼并、重组。

（四） 有效发挥WTO对于重庆汽车产业的积极作用

重庆在面对WTO的机遇与挑战时，要利用好WTO有关规则，对重庆汽车产业进行保护和扶持。借鉴发达国家的成功经验，按照WTO有关条款及成员国(方)的惯例做法，对本地汽车企业制定相应政策，实施有效的保护和扶持。随着进口关税进一步下调，配额不断增加，国内汽车贸易体系对外国公司的进一步开放，汽车企业将面临着更加严峻的竞争形势。提高自主创新能力，加强产品核心竞争力，加大与国外大型知名汽车企业的合作。在新形势与新挑战下，创造出适合重庆汽车产业发展的新道路。

七、结束语

在未来的5到10年间，重庆汽车的发展目标应是走新型工业化道路，提高核心竞争力，力争把重庆建设成为我国西部“汽车名城”。为了达到这一目标必须坚定不移地以改革为动力，以龙头企业为带动力量，以企业自主创新为依托，全面加快汽车产业的发展。重庆汽车要真正崛起，必须引进人才、整合产业链、走研发创新之路。依托重庆钢铁、綦江齿轮等优势产业，重庆从零部件，到汽车整车生产企业的互相渗透和整合，才是必然之路。

参考文献

[1]孙云利，万建平，黄光辉．转变点分析在我国外贸安全中的应用．统计创新与实践．湖北省统计学会．武汉：崇文书局，2003：4l一48，

[2]孙云利，万建平，黄光辉．转变点识别的稳定性．大学数学，2004(2)．(拟发表)

[3]孙云利，万建平，黄光辉．非参数回归分析方法及其在经济模型中的应用．应用数学，2003，16(增)：176．179．

[4]杨虎，刘琼荪，钟波 . 数理统计.北京： 高等教育出版社，2004.10

[5]王学民.应用多元分析.上海：上海财经大学出版社，2009.8

[6] 董建华．汽车行业：2006年度基础分析报告及建议[EB／OL]．2007.5

[7] 毕建平，王新莉，邵风翔，马建辉．中国民用汽车拥有量时间序列动态模型研究[J]．河南科学，2004(12)：847—849．

[8] 郑祖康．假设检验的稳定性[J]．中国科学(A辑)，1997，27(9)：782—792．

[9] 王彦森. 中国汽车产业集群化发展的可行性分研究[J].技术经济，2004.7:6-9

[10] 熊文才.对重庆汽车工业发展的思考[J] .重庆工学院学报 ,2002(8) :85 - 90.

[11] 郝云宏.中国汽车工业国际竞争力分析[J] .商业经济与管理 ,2003(3) :4 - 9.

[12] 徐小钦,谭利华. 汽车产业提升自主创新能力研究——以重庆为例. 重庆大学学报 [J].2007年第13卷第2期

[13] 高　韵,罗有贤. 重庆汽车产业发展战略[J].重庆工学院学报, 2008.6

[14] 詹　乔,陈静宇.加入WTO对重庆汽车产业的影响及对策[J] .重庆大学学报, 2000年第6卷第2期

[15] 重庆市政府.重庆市人民政府办公厅关于转发重庆市2001年城市道路交通畅通工程实施意见的通知(渝办发[2001]23).http://www.cq.gov.cn/zwgk/zfgw/51971.htm. 2004年06月11日

[16] 重庆市政府.重庆市人民政府办公厅转发市公安局等部门关于重庆市2002年城市道路畅通工程实施意见补充方案的通知.http://www.cq.gov.cn/zfgbyxwfb/zfgb/gb2002/15qi/57003.htm .2004年05月14日

[17] 重庆市政府.重庆市人民政府关于认真做好2002年市重点建设项目和重大前期工作项目有关工作的通知. http://www.cq.gov.cn/zfgbyxwfb/zfgb/gb2002/dijiuqi/56724.htm. 2004年05月14日

附录1

一、 回归模型、主成分分析、主成分回归模型程序

car=read.csv("car.csv") %（在R软件中读取改变成csv文件的Excel中的数据）

attach(car)

opar<-par(mfrow=c(2,2), mar=c(5,4,1,1))

plot(x1, y,type='o'); plot(x2, y,type='o'); plot(x3, y,type='o'); plot(x4, y,type='o') par(opar)> lm.sol<-lm(y~x1+x2+x3+x4, data=car)

summary(lm.sol)

结果如下：

Call:

lm(formula = y ~ x1 + x2 + x3 + x4, data = car)

Residuals:

Min 1Q Median 3Q Max

-0.27948 -0.11034 -0.02974 0.11112 0.40336

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 188.65723 182.24760 1.035 0.325

x1 -0.13805 0.58320 -0.237 0.818

x2 -0.03798 6.93596 -0.005 0.996

x3 - 5.87879 5.59459 1.051 0.318

x4 -57.95460 55.03946 - 1.053 0.317

Residual standard error: 0.2187 on 10 degrees of freedom

Multiple R-squared: 0.9316, Adjusted R-squared: 0.9043

F-statistic: 34.06 on 4 and 10 DF, p-value: 8.46e-06

> car.pr<-princomp(~x1+x2+x3+x4, data=car, cor=T)

> summary(car.pr, loadings=TRUE)

结果：

Importance of components:

Comp.1 Comp.2 Comp.3 Comp.4

Standard deviation 1.968206 0.34596826 0.073701602 0.0322208910

Proportion of Variance 0.968459 0.02992351 0.001357982 0.0002595465

Cumulative Proportion 0.968459 0.99838247 0.999740454 1.0000000000

Loadings:

Comp.1 Comp.2 Comp.3 Comp.4

x1 0.484 0.874

x2 0.506 -0.239 -0.449 0.697

x3 0.505 -0.294 -0.384 -0.715

x4 0.504 -0.305 0.807

> pre<-predict(car.pr)

> car$Z<-pre[,1]

> lm.sol<-lm(y~Z, data=car)

> summary(lm.sol)

Call:

lm(formula = y ~ Z, data = car)

Residuals:

Min 1Q Median 3Q Max

-0.19910 -0.13114 -0.07941 0.04380 0.52884

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 5.27851 0.05892 89.58 < 2e-16 ***

Z 0.32973 0.02994 11.01 5.8e-08 ***

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2282 on 13 degrees of freedom

Multiple R-squared: 0.9032, Adjusted R-squared: 0.8958

F-statistic: 121.3 on 1 and 13 DF, p-value: 5.802e-08

> beta<-coef(lm.sol); A<-loadings(car.pr)

> x.bar<-car.pr$center; x.sd<-car.pr$scale

> coef<-beta[2]*A[,1]/x.sd

> beta0 <- beta[1]- sum(x.bar * coef)

> c(beta0, coef)

(Intercept) x1 x2 x3 x4

-54.6081213 0.5962799 0.9116767 0.7234347 14.6212501

用MATLAB软件计算转变点

1、计算SIC(15)

X=[

1 4.353262 3.641021 3.050403 3.477377;1 4.429623 3.70096 3.118965 3.480405;

1 4.432087 3.724444 3.178905 3.483291; 1 4.434729 3.73582 3.204766 3.485677;

1 4.44849 3.765552 3.220944 3.487469;1 4.482216 3.790728 3.253096 3.490112;

1 4.486487 3.817717 3.295976 3.491069;1 4.492201 3.859623 3.348861 3.493295;

1 4.497026 3.908145 3.407513 3.495558;1 4.509794 3.9647763.482099 3.497514;

1 4.992191 4.010469 3.540044 3.500944;15.001297 4.063324 3.591869 3.504997;

1 5.019967 4.137204 3.669887 3.509917;1 5.0359584.196141 3.762953 3.512824; 1 5.045131 4.235304 3.814914 3.515292]；

Y=[4.357649 4.424163 4.541217 4.63794 4.663739 4.762198 4.91598 5.015804 5.739714 5.793644 5.93406 6.01181 6.051865 6.107711 6.220127]'；

b=inv(X'*X)*X'*Y；

SIC(15)=15*log[(Y-X*b)'*(Y-X*b)]+15*((log2*pi)+1)-9*log(15)

ans=

7.1331

算个各转变点的程序：

K=5

X3=[1 4.353262 3.641021 3.050403 3.477377;1 4.429623 3.70096 3.118965 3.480405;1 4.432087 3.724444 3.178905 3.483291; 1 4.434729 3.735826 3.204766 3.485677;1 4.44849 3.765552 3.220944 3.487469]；

X4=[1 4.482216 3.790728 3.253096 3.490112;1 4.486487 3.817717 3.295976 3.491069;1 4.492201 3.859623 3.348861 3.493295;1 4.497026 3.908145 3.407513 3.495558;1 4.509794 3.964776 3.482099 3.497514;1 4.992191 4.010469 3.540044 3.500944;1 5.001297 4.063324 3.591869 3.504997; 1 5.019967 4.137204 3.669887 3.509917;1 5.035958 4.196141 3.762953 3.512824;1 5.045131 4.235304 3.814914 3.515292]；

Y3=[4.357649 4.424163 4.541217 4.63794 4.663739]'；

Y4=[4.762198 4.91598 5.015804 5.739714 5.793644 5.93406 6.01181 6.051865 6.107711 6.220127]'；

b3=inv(X3'*X3)*X3'*Y3；

b4=inv(X4'*X4)*X4'*Y4；

SIC(5)=15*log[(Y3-X3*b3)'*(Y3-X3*b3)+(Y4-X4*b4)'*(Y4-X4*b4)]+15*(log(2*pi)+1)-4*log(15)

SIC =

0 0 0 0 2.6251 0 0 0 0 0 0 0 0 0 7.1331

K=6

X5=[1 4.353262 3.641021 3.050403 3.477377;1 4.429623 3.70096 3.118965 3.480405;1 4.432087 3.724444 3.178905 3.483291;1 4.434729 3.735826 3.204766 3.485677;1 4.44849 3.765552 3.220944 3.487469;1 4.482216 3.790728 3.253096 3.490112];

X6=[

1 4.486487 3.817717 3.295976 3.491069;1 4.492201 3.859623 3.3488613.493295;

1 4.497026 3.908145 3.407513 3.495558;1 4.509794 3.964776 3.4820993.497514;

1 4.992191 4.010469 3.540044 3.500944;1 5.001297 4.063324 3.5918693.504997;

1 5.019967 4.137204 3.669887 3.509917;1 5.035958 4.1961413.762953 3.512824;

1 5.045131 4.235304 3.814914 3.515292];

Y5=[4.357649 4.424163 4.541217 4.63794 4.663739 4.762198]';

Y6=[4.91598 5.015804 5.739714 5.793644 5.93406 6.01181 6.051865 6.107711 6.220127]';

b5=inv(X5'*X5)*X5'*Y5;

b6=inv(X6'*X6)*X6'*Y6;

SIC(6)=15*log[(Y5-X5*b5)'*(Y5-X5*b5)+(Y6-X6*b6)'*(Y6-X6*b6)]+15*(log(2*pi)+1)-4*log(15)

SIC =

0 0 0 0 2.6251 2.2013 0 0 0 0 0 0 0 0 7.1331

K=7

X9=[1 4.353262 3.641021 3.050403 3.477377;1 4.429623 3.70096 3.118965 3.480405;1 4.432087 3.724444 3.178905 3.483291;1 4.434729 3.735826 3.204766 3.485677;1 4.44849 3.765552 3.220944 3.487469;1 4.482216 3.790728 3.253096 3.490112; 1 4.486487 3.817717 3.295976 3.491069];

X10=[ 1 4.492201 3.859623 3.348861 3.493295;1 4.497026 3.908145 3.407513 3.495558; 1 4.509794 3.964776 3.482099 3.497514;1 4.992191 4.010469 3.540044 3.500944;1 5.001297 4.063324 3.591869 3.504997;1 5.019967 4.137204 3.669887 3.509917;1 5.035958 4.196141 3.762953 3.512824;1 5.045131 4.235304 3.814914 3.515292];

Y9=[4.357649 4.424163 4.541217 4.63794 4.663739 4.762198 4.91598]';

Y10=[5.015804 5.739714 5.793644 5.93406 6.01181 6.051865 6.107711 6.220127]';

b9=inv(X9'*X9)*X9'*Y9;

b10=inv(X10'*X10)*X10'*Y10;

SIC(7)=15*log[(Y9-X9*b9)'*(Y9-X9*b9)+(Y10-X10*b10)'*(Y10-X10*b10)

]+15*(log(2*pi)+1)-4*log(15)

SIC =

0 0 0 0 2.6251 2.2013 1.6765 0 0 0 0 0 0 0 7.133

K=8

X1=[1 4.353262 3.641021 3.050403 3.477377;1 4.429623 3.70096 3.118963.480405;1 4.432087 3.724444 3.178905 3.483291;1 4.434729 3.735826 3.204766 3.485677;1 4.44849 3.765552 3.220944 3.487469;1 4.482216 3.790728 3.253096 3.490112; 1 4.486487 3.817717 3.295976 3.491069;1 4.492201 3.8596233.3488613.493295];

Y1=[4.357649 4.424163 4.541217 4.63794 4.663739 4.762198 4.91598 5.015804]';

b1=inv(X1'*X1)*X1'*Y1;

X2=[1 4.497026 3.908145 3.407513 3.495558; 1 4.509794 3.964776 3.482099 3.497514;1 4.992191 4.010469 3.540044 3.500944;1 5.001297 4.063324 3.591869 3.504997;1 5.019967 4.137204 3.669887 3.509917;1 5.035958 4.196141 3.762953 3.512824;1 5.045131 4.235304 3.814914 3.515292];

Y2=[5.739714 5.793644 5.93406 6.01181 6.051865 6.107711 6.220127]';

b2=inv(X2'*X2)*X2'*Y2;

SIC(8)=15*log[(Y1-X1*b1)'*(Y1-X1*b1)+ Y2-X2*b2)'*(Y2-X2*b2)]+15*((log2*pi)+1)-4*log(15)

SIC (8)=

0 0 0 0 2.6251 2.2013 1.6765 -40.5024 0 0 0 0 0 0 7.133

K=9

X9=[1 4.353262 3.641021 3.050403 3.477377;1 4.429623 3.70096 3.118965 3.480405;1 4.432087 3.724444 3.178905 3.483291;1 4.434729 3.735826 3.204766 3.485677;1 4.44849 3.765552 3.220944 3.4874691;1 4.482216 3.790728 3.253096 3.490112; 1 4.486487 3.817717 3.295976 3.491069;1 4.492201 3.85962 3.348861 3.493295;1 4.497026 3.908145 3.407513 3.495558];

X10=[1 4.509794 3.964776 3.482099 3.497514;1 4.992191 4.010469 3.540044 3.500944;1 5.001297 4.063324 3.591869 3.504997;1 5.019967 4.137204 3.669887 3.509917; 1 5.035958 4.196141 3.762953 3.512824;1 5.045131 4.235304 3.814914 3.515292];

Y9=[4.357649 4.424163 4.541217 4.63794 4.663739 4.762198 4.91598 5.015804 5.739714]';

Y10=[5.793644 5.93406 6.01181 6.051865 6.107711 6.220127]';

b9=inv(X9'*X9)*X9'*Y9;

b10=inv(X10'*X10)*X10'*Y10;

SIC(9)=15*log[(Y9-X9*b9)'*(Y9-X9*b9)+(Y10-X10*b10)'*(Y10-X10*b10)]

+15*(log(2*pi)+1)-4*log(15)

SIC =

0 0 0 0 2.6251 2.2013 1.6765 -40.5024 -7.2587 0 0 0 0 0 7.1331

K=10

X13=[1 4.353262 3.641021 3.050403 3.477377;1 4.429623 3.70096 3.118965 3.480405;1 4.432087 3.724444 3.178905 3.483291; 1 4.434729 3.735826 3.204766 3.485677;1 4.44849 3.765552 3.220944 3.4874691;1 4.482216 3.790728 3.253096 3.490112; 1 4.486487 3.817717 3.295976 3.491069;1 4.492201 3.859623 3.348861 3.493295;1 4.497026 3.908145 3.407513 3.495558; 1 4.509794 3.964776 3.482099 3.497514];

X14=[ 1 4.99219 1 4.010469 3.540044 3.500944;1 5.001297 4.063324 3.591869 3.504997;1 5.019967 4.137204 3.669887 3.509917;1 5.035958 4.1961413.762953 3.512824;1 5.045131 4.235304 3.814914 3.515292];

Y13=[4.357649 4.424163 4.541217 4.63794 4.663739 4.762198 4.91598 5.015804 5.739714 5.793644]';

Y14=[5.93406 6.01181 6.051865 6.107711 6.220127]';

b13=inv(X13'*X13)*X13'*Y13;

b14=inv(X14'*X14)*X14'*Y14;

SIC(10)=15*log[(Y13-X13*b13)'*(Y13-X13*b13)+(Y14-X14*b14)'*(Y14-X14*b14)]

+15*(log(2*pi)+1)-4*log(15)

SIC =

0 0 0 0 2.6251 2.2013 1.6765 -40.5024 -7.2587 -3.8270 129.2805 0 0 0 7.1331

K=11

X11=[1 4.353262 3.641021 3.050403 3.477377;1 4.429623 3.70096 3.118965 3.480405;1 4.432087 3.724444 3.178905 3.483291; 1 4.434729 3.735826 3.204766 3.485677;1 4.44849 3.765552 3.220944 3.487469;1 4.482216 3.790728 3.253096 3.490112;1 4.486487 3.817717 3.295976 3.491069;1 4.492201 3.859623 3.348861 3.493295;1 4.497026 3.908145 3.407513 3.495558;1 4.509794 3.964776 3.482099 3.497514;1 4.992191 4.010469 3.540044 3.500944];

X12=[1 5.001297 4.063324 3.591869 3.504997; 1 5.019967 4.137204 3.669887 3.509917;1 5.035958 4.196141 3.762953 3.512824;1 5.045131 4.235304 3.814914 3.515292];

Y11=[4.357649 4.424163 4.541217 4.63794 4.663739 4.762198 4.91598 5.015804 5.739714 5.793644 5.93406]';

Y12=[6.01181 6.051865 6.107711 6.220127]';

b11=inv(X11'*X11)*X11'*Y11;

b12=inv(X12'*X12)*X12'*Y12;

SIC(11)=15*log[ (Y11-X11*b11)'*(Y11-X11*b11)+(Y12-X12*b12)'*(Y12-X12*b12)]+15*(log(2*pi)+1)-4*log(15)

SIC =

0 0 0 0 2.6251 2.2013 1.6765 -40.5024 -7.2587 -3.8270 129.2805 0 0 0 7.1331

得出转变点后进行的回归分析程序

1、转变前：

car=read.csv("car1.csv")

attach(car)

opar<-par(mfrow=c(2,2), mar=c(5,4,1,1))

plot(x1, y,type='o'); plot(x2, y,type='o'); plot(x3, y,type='o'); plot(x4, y,type='o')

par(opar)

m.sol<-lm(y~x1+x2+x3+x4, data=car)

summary(lm.sol)

car.pr<-princomp(~x1+x2+x3+x4, data=car, cor=T)

summary(car.pr, loadings=TRUE)

pre<-predict(car.pr)

car$Z<-pre[,1]

lm.sol<-lm(y~Z, data=car)

summary(lm.sol)

beta<-coef(lm.sol); A<-loadings(car.pr)

x.bar<-car.pr$center; x.sd<-car.pr$scale

coef<-beta[2]*A[,1]/x.sd

beta <- beta[1]- sum(x.bar * coef)

c(beta, coef)

转变后：

car=read.csv("car2.csv")

attach(car)

opar<-par(mfrow=c(2,2), mar=c(5,4,1,1))

plot(x1, y,type='o'); plot(x2, y,type='o'); plot(x3, y,type='o'); plot(x4, y,type='o')

par(opar)

m.sol<-lm(y~x1+x2+x3+x4, data=car)

summary(lm.sol)

car.pr<-princomp(~x1+x2+x3+x4, data=car, cor=T)

summary(car.pr, loadings=TRUE)

pre<-predict(car.pr)

car$Z<-pre[,1]

lm.sol<-lm(y~Z, data=car)

summary(lm.sol)

beta<-coef(lm.sol); A<-loadings(car.pr)

x.bar<-car.pr$center; x.sd<-car.pr$scale

coef<-beta[2]*A[,1]/x.sd

beta <- beta[1]- sum(x.bar * coef)

c(beta, coef)

本文来源：https://www.2haoxitong.net/k/doc/dcde6262af1ffc4ffe47acea.html