epoh老师，利用贝叶斯估计上述模型在R的那个packages
可以实现，我只找到了包含threshold AR的，没有找到
不包含threshold的AR模型。

try package 'MCMCpack'
function MCMCregress()
  Markov Chain Monte Carlo for Gaussian Linear Regression

这个package我看过了，大致模型应该是一致的，当时考虑上述包含滞后性的AR模型
是不是与Gaussian Linear Regression完全一致。
谢谢epoh老师了，我再看下这个package。

# generate an AR(1) with mean 50
set.seed(66)      # so you can reproduce these results
x = arima.sim(list(order=c(1,0,0), ar=.9), n=100) + 50   
mean(x)  #  [1] 50.60668    the sample mean is close

y=x[2:100]
y1=x[1:99]

#######
lm=lm(y~y1)
summary(lm)

#Coefficients:
#            Estimate Std. Error t value Pr(>|t|)   
#(Intercept)  4.80587    2.18935   2.195   0.0305 *  
#y1              0.90523    0.04323  20.941   <2e-16 ***

#########
library(MCMCpack)
posterior <- MCMCregress(y~y1, verbose=1000)
plot(posterior)
raftery.diag(posterior)
summary(posterior)

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

              Mean      SD  Naive SE Time-series SE
(Intercept) 4.8182 2.21104 0.0221104       0.023461
y1             0.9050 0.04368 0.0004368       0.000467
sigma2      0.9164 0.13541 0.0013541       0.001512

谢谢epoh老师，明白了。

各位前辈，我想问下，为什么我估计出的LSTAR模型中，t值显示为NA?
> mod=lstar(y,m=4,d=1,thVar=x)
Using maximum autoregressive order for low regime: mL = 4
Using maximum autoregressive order for high regime: mH = 4
Using only first 77 elements of thVar
Performing grid search for starting values...
Starting values fixed: gamma =  100 , th =  0.0319328 ; SSE =  0.0642117
Grid search selected lower/upper bound gamma (default [1,100]).
                                          Might try to widen bound with arg: 'starting.control=list(gammaInt=c(1,200))'
Optimization algorithm converged
Optimized values fixed for regime 2  : gamma =  100.001 , th =  0.1487109 ; SSE =  0.05262752
> mod

Non linear autoregressive model

LSTAR model
Coefficients:
Low regime:
    const1     phi1.1     phi1.2     phi1.3     phi1.4
0.0348664 -0.1229966 -0.1646415 -0.1743049  0.5060107

High regime:
     const2      phi2.1      phi2.2      phi2.3      phi2.4
   666.6409 -10646.2203  -1317.4663   8211.4405   4816.7396

Smoothing parameter: gamma = 100

Threshold
Variable: external
Value: 0.1487
> summary(mod)

Non linear autoregressive model

LSTAR model
Coefficients:
Low regime:
    const1     phi1.1     phi1.2     phi1.3     phi1.4
0.0348664 -0.1229966 -0.1646415 -0.1743049  0.5060107

High regime:
     const2      phi2.1      phi2.2      phi2.3      phi2.4
   666.6409 -10646.2203  -1317.4663   8211.4405   4816.7396

Smoothing parameter: gamma = 100

Threshold
Variable: external
Value: 0.1487

Residuals:
        Min          1Q      Median          3Q         Max
-0.06801287 -0.01534197  0.00094419  0.01345869  0.06109822

Fit:
residuals variance = 0.0006497,  AIC = -570, MAPE = 110.3%

Coefficient(s):
          Estimate  Std. Error  t value Pr(>|z|)
const1  3.4866e-02          NA       NA       NA
phi1.1 -1.2300e-01          NA       NA       NA
phi1.2 -1.6464e-01          NA       NA       NA
phi1.3 -1.7430e-01          NA       NA       NA
phi1.4  5.0601e-01          NA       NA       NA
const2  6.6664e+02          NA       NA       NA
phi2.1 -1.0646e+04          NA       NA       NA
phi2.2 -1.3175e+03          NA       NA       NA
phi2.3  8.2114e+03          NA       NA       NA
phi2.4  4.8167e+03          NA       NA       NA
gamma   1.0000e+02          NA       NA       NA
th      1.4871e-01          NA       NA       NA

Non-linearity test of full-order LSTAR model against full-order AR model
F = 3.8863 ; p-value = 0.0066813

Threshold
Variable: external警告信息：
In summary.lstar(mod) : singular Hessian

>

epoh老师您好，好久不见了！我有个问题请教您，之前您教我做的那个脉冲影响分析，为什么我给STAR模型之后做出来相同大小的正的脉冲和负的脉冲值是大小相等，方向相反的呢？而且1个单位的脉冲和2个单位的脉冲也是2倍关系？一个外文文献中讲到，非线性模型中的的脉冲是不具有线性和对称性的？可是为什么我做出来的是有线性和对称性的呢？难道是模型出错了？期待您的指教，谢谢啦！

请先参考81楼我给南冰兄的asym.jpeg
建议你朝两个方向探讨:

1.南冰兄的quartc.txt, "qcpi"
   与你的x11.txt, "y"数据大小,有些差距.

2.改变一下 x11.txt, "s"试试

老师您好！我看了81楼您发的程序，和之前您发给我的那个脉冲的比较了一下，好像没发现怎么可以处理我的问题的办法，我又做了一遍，还是对称的。那您说我之前做的LSTAR模型的脉冲，同样是两个标准差但是符号相反的脉冲函数值会是大小相等，方向相反，呈对称嘛？应该不会是模型出问题了吧？不好意思，我的问题太多了，给您添麻烦了！

我的意思是你的数据"y",

处理过程有没笔误,因为看不出有Business Cycle.

底下我提供Terasvirta JASA 1994

Terasvirta(1994), "Specification, Estimation and Evaluation of

                            Smooth Transition Autoregressive Models", JASA, vol 89,

data "lynx",及 "ipgrowth"的图型供你参考.

#####

library(tsDyn)
#data lynx

#data ipgrowth
ip =read.table("ipgrowth.txt", header = TRUE);
ipgrowth =ip$ipgrowth

#your data y
new=read.table("x11.txt",header=TRUE)
y=new$y

#plot
par(mfrow = c(2, 2))
plot(log10(lynx))
lines(log10(lynx),col='red')

plot(ipgrowth)
lines(ipgrowth,col='red')

plot(y)
lines(y,col='red')

哦，我明白您的意思了！我对Y是进行了零均值化和X-12季节调整的，可能是期中的这种Cycle就显示不出来了，我回头再试试！还有个问题，对于下面做出来的脉冲响应函数，是不是可以说，受到外来1个单位的冲击时，Y在将来一期时波动是0.4，然后随着时间不断减小呢？就是这个脉冲响应是不是可以作为Y的一个预测呢？谢谢您了epoh老师！

1.脉冲响应函数的解释可以参考
  易丹辉2008数据分析与EViews应用.pdf
  chap 8,page 227/408
  https://bbs.pinggu.org/thread-998584-1-1.html

2.Y的一个预测,可以使用 function predict()
  ## S3 method for class 'nlar':
     predict(object, newdata, n.ahead, simulate=FALSE, ...)
  请参考Nonlinear autoregressive time series models in R.pdf
  page 21/35,4.3 Out-of-sample forecasting
  http://cran.r-project.org/web/packages/tsDyn/index.html

好的，老师我明白了，谢谢您了啊！

epoh老师，您好！好久不见啊，想请教您个问题，就是在STAR模型中，我看好些外文文献都把协整引进来了，多变量引入协整的时候，是把协整的关系式放在转换变量的位置上呢还是放在线性部分呢？我有点疑惑，没看多懂，那这种情况用R软件能实现估计吗？谢谢您了！

这跟一般的vecm model,观念接近.
请参考文献:货币—产出lstvecm.pdf  
page 3/13 公式(1)及对照 page7/13 Table 3的各项系数,就会明白

目前这个模型,似乎没有现成的R package可供使用.

好的，谢谢epoh老师！我知道了，一定好好阅读这个文章！

Oil Prices and the Business Cycle.pdf
  

请问在gauss编程中，我想调用MSVARlib-v2.0中的库函数(.lib)和source函数（.src),
请问这些.lib和.src函数如何调用？有什么区别和联系？
另外，我发现一些gauss程序中，在使用了.lib库文件后，还用.dll文件，这两个文件有什么区别和联系？
如何使用？

请问epoh老师，
请问在gauss编程中，我想调用MSVARlib-v2.0中的库函数(.lib)和source函数（.src),
请问这些.lib和.src函数如何调用？有什么区别和联系？
另外，我发现一些gauss程序中，在使用了.lib库文件后，还用.dll文件，这两个文件有什么区别和联系？
如何使用？
能给个具体 的例子吗？
非常感谢！

抱歉!zhangtao兄,最近较忙,晚回了.
刚抽空看了下
MSVARlib-v2.0并没有你说的库函数(.lib)和source函数（.src),
作者只用.prg,简单明了
要用.prg只要
  #include XXX.prg;
要用.src只要
  #include XXX.src;
范例:
C:\gauss\MSVAR\MSlib\Examples\Model1.prg
   #include  C:\Gauss\MSVAR\MSLIB\MSVAR_Call.prg;  
    .......
    .......

C:\gauss\MSVAR\MSlib\MSVAR_Call.prg;  再呼叫底下函数
/*==========================================================================================*/
/* MSVAR_Call : calls different subroutines need to run the main program                    */
/*==========================================================================================*/

/*==============Loading, setting and statistics files================*/

#include C:\Gauss\MSVAR\MSlib\MSVAR_load.prg;
#include C:\Gauss\MSVAR\MSlib\MSVAR_Setdatation.prg
#include C:\Gauss\MSVAR\MSlib\MSVAR_Setsample.prg;
#include C:\Gauss\MSVAR\MSlib\MSVAR_moment.prg;


/*============Intialisation and Transformation-constraint files ================*/

#include C:\Gauss\MSVAR\MSlib\MSVAR_Initsort.prg;
#include C:\Gauss\MSVAR\MSlib\MSVAR_Init.prg;   
#include C:\Gauss\MSVAR\MSlib\MSVAR_VecMat.prg;


/*============Filters and smoothers files ================*/

#include C:\Gauss\MSVAR\MSlib\MSVAR_FiltHmm.prg;
#include C:\Gauss\MSVAR\MSlib\MSVAR_SmoHmm.prg;

/*============Maximum likelihood files================*/

#include C:\Gauss\MSVAR\MSlib\MSVAR_MaxHmm.prg;
#include C:\Gauss\MSVAR\MSlib\MSVAR_stderr.prg;

/* ============ Export files ==========================*/

#include C:\Gauss\MSVAR\MSlib\MSVAR_Output.prg;  


/****************/
gauss 和R一样,可以调用C,Fortran,..等函数,所以会有Dynamic Link Library
zhangtao兄有兴趣可以参考底下网页:
  A GAUSS/DLL Implementation of Alan Genz’s Fortran Packages
  for Computing Multivariate t CDF
  http://www.thierry-roncalli.com/download/gauss-mvt.pdf
下载gauss-mvt.zip参考

epoh老师，您好！
         为什么附件 中的examples文件我运行总是提示以下错误,为什么？如何更正？
非常感谢！
          另外，我还是没有明白.lib文件和.dll文件有什么区别和联系。
          另外，编程中（例如java）中的多态是什么？我总是不明白。
Error in ==> Example_MS_Regress_Fit_using_constCoeff at 35
[Spec_Out]=MS_Regress_Fit(dep,indep,k,S,advOpt); % Estimating the model
??? Error using ==> m_Files\checkSize_constCoeff.m
The argument constcoeff.nS_Param{iEq} should be a vector with number of rows equal to the number of non switching parameters (indep matrix)
Error in ==> m_Files\checkSize_constCoeff.m at 24
            error('The argument constcoeff.nS_Param{iEq} should be a vector with number of rows equal to the number of non switching parameters (indep matrix)')
Error in ==> m_Files\MS_Regress_Fit.m at 39
checkSize_constCoeff();
Error in ==> Example_MS_Regress_For at 20
[Spec_Out]=MS_Regress_Fit(dep,indep,k,S,advOpt); % Estimating the model
??? Error using ==> fmincon
FMINCON cannot continue because user supplied objective function failed with the following error:
Variable 'indep_nS' is used as a command function.
Error in ==> m_Files\MS_Regress_Fit.m at 85
[param]=fmincon(@(param)MS_Regress_Lik(dep,indep_nS,indep_S,param,k,S,advOpt,dispOut),param0, ...
Error in ==> Example_MS_Regress_Simul_and_Fit_2_States at 54
[Spec_Output]=MS_Regress_Fit(dep,indep,k,S,advOpt);
??? Error using ==> fmincon
FMINCON cannot continue because user supplied objective function failed with the following error:
Variable 'indep_nS' is used as a command function.
Error in ==> m_Files\MS_Regress_Fit.m at 85
[param]=fmincon(@(param)MS_Regress_Lik(dep,indep_nS,indep_S,param,k,S,advOpt,dispOut),param0, ...
Error in ==> Example_MS_Regress_Simul_and_Fit_3_States at 56
[Spec_Output]=MS_Regress_Fit(dep,indep,k,S,advOpt); % fit the model
??? Error using ==> fmincon
FMINCON cannot continue because user supplied objective function failed with the following error:
Variable 'indep_nS' is used as a command function.
Error in ==> m_Files\MS_Regress_Fit.m at 85
[param]=fmincon(@(param)MS_Regress_Lik(dep,indep_nS,indep_S,param,k,S,advOpt,dispOut),param0, ...
Error in ==> Example_MS_Regress_Simul_Fit_2_States_Multivar at 65
[Spec_Out]=MS_Regress_Fit(dep,indep,k,S,advOpt); % Estimating the model
>>

c:\gauss10\src\mvt.src(44) : error G0432 : 'c:\gauss10\dlib\mvt1.dll' : Error opening dynamic library
LoadLibrary error 126 occurred

c:\gauss10\src\mvt.src(44) : error G0432 : 'c:\gauss10\dlib\mvt2.dll' : Error opening dynamic library
为什么附件中的gauss程序在gauss10中运行会出现以错误？

1.程序可以正常运转
2."java的多态",不太懂这句话的意思,请用英文表达
3.抱歉!目前我的computer没安装Gauss,
   Error opening dynamic library
    可能是你没有把 mvt1.dll,mvt2.dll
    放进dlib folder
%%%%%%%%%%
set path\add with subfolders\MS_Regress_FEX
>> Example_MS_Regress_Fit_using_constCoeff
Calculating Standard Error Vector...


***** Numerical Optimization Converged *****

Final log Likelihood: 2425.7305
Number of estimated parameters: 6
Type of Switching Model: Univariate
Distribution Assumption -> Normal
Method SE calculation -> 1

***** Final Parameters for Equation #1 *****

---> Non Switching Parameters <---

Non Switching Parameter for Equation #1, Indep column 1
     Value:                0.0000
     Std Error (p. value): Inf (1.00)

--->   Switching Parameters (Distribution Parameters)  <---

State 1
    Model's Variance:      0.000100
    Std Error (p. value):  NaN (NaN)
State 2
    Model's Variance:      0.000459
    Std Error (p. value):  0.0000 (0.00)

--->   Switching Parameters (Regressors)  <---

Switching Parameters for Equation #1 - Indep column 2

State 1
   Value:                0.0000
   Std Error (p. value): NaN (NaN)
State 2
   Value:                0.5339
   Std Error (p. value): 0.0267 (0.00)

Switching Parameters for Equation #1 - Indep column 3

State 1
   Value:                0.9869
   Std Error (p. value): 0.0898 (0.00)
State 2
   Value:                0.0000
   Std Error (p. value): NaN (NaN)

---> Transition Probabilities Matrix (std. error, p-value) <---

      0.95 ( NaN, NaN)   0.01 (0.00,0.22)   
      0.05 ( NaN, NaN)   0.99 (0.04,0.00)   

---> Expected Duration of Regimes <---

     Expected duration of Regime #1: 20.00 time periods
     Expected duration of Regime #2: 181.35 time periods

我刚重新安装GAUSS9.0,
试了可以运行出结果:
     Upper    cdfmvt      cdft     diff.
=========================================
    -3.000     0.048     0.048     0.000
    -2.500     0.065     0.065     0.000
    -2.000     0.092     0.092     0.000
    -1.500     0.136     0.136     0.000
    -1.000     0.211     0.211     0.000
    -0.500     0.333     0.333     0.000
     0.000     0.500     0.500     0.000
     0.500     0.667     0.667     0.000
     1.000     0.789     0.789     0.000
     1.500     0.864     0.864     0.000
     2.000     0.908     0.908     0.000
     2.500     0.935     0.935     0.000
     3.000     0.952     0.952     0.000


     Upper  cdfmvtnc    cdftnc     diff.
=========================================
    -3.000     0.002     0.002     0.000
    -2.500     0.003     0.003     0.000
    -2.000     0.005     0.005     0.000
    -1.500     0.008     0.008     0.000
    -1.000     0.014     0.014     0.000
    -0.500     0.029     0.029     0.000
     0.000     0.067     0.067     0.000
     0.500     0.152     0.152     0.000
     1.000     0.287     0.287     0.000
     1.500     0.436     0.436     0.000
     2.000     0.566     0.566     0.000
     2.500     0.666     0.666     0.000
     3.000     0.740     0.740     0.000
请务必依照第一页指示,放置文件

target path          readme.mvt
target path\dlib     Dynamic link libraries
target path\lib      Library file
target path\mvt      Examples and tutorial files
target path\mvt\dlls Fortran/C source code files
target path\src      GAUSS source code files

你的target path就是
c:\gauss10

质量很高的讨论贴 强烈赞一个

今天恰巧有人问我底下这篇文献,
cyclical component如何画出来.
Markov机制转换的状态空间模型及其在我国经济周期中的应用研
  http://wenku.baidu.com/view/0cdd2961a45177232f60a2fb.html

我想到这篇文献对你应有帮助
Nonlinearities over the business cycle_banco
       http://www.bportugal.pt/en-US/BdP%20Publications%20Research/WP200309.pdf
这篇文献就是用smooth transition autoregressive (STAR) model

太感谢您了！好的，我下载了一定好好看，谢谢啦

epoh老师确实懂得很多，不知道能给出stata的estar程序就好了！

good

你好：我用你的数据和方法，得到了结果确不同与你的，我的tsDyn是9.o版的。结果如下：
>svpdx <- read.table("data.txt", header = TRUE)
>svpdx
    cpi  rpi          m2
1   3.4  2.9  0.26524514
2   6.4  5.4  0.31278198
3  14.7 13.2  0.37310154
4  24.1 21.7  0.34529154
5  17.1 14.8  0.29467111
6   8.3  6.1 -0.03867458
7   2.8  0.8  0.34668927
8  -0.8 -2.6  0.17245702
9  -1.4 -3.0  0.14411363
10  0.4 -1.5 -0.31433668
11  0.7 -0.8  0.87572994
12 -0.8 -1.3  0.16337324
13  1.2 -0.1  0.19679602
14  3.9  2.8  0.18161695
15  1.8  0.8  0.14387323
16  1.5  1.0  0.18886459
17  4.8  3.8  0.15787746
18  5.9  5.9  0.18875571
19 -0.7 -1.2  0.18736354
> x=svpdx$cpi
> y=svpdx$rpi
> z=svpdx$m2
> library(tsDyn)
>ndx.lstar <- lstar(x, m=3,d=1, thVar=z,control=list(maxit=3000))
（使用此命令得到的结果与你的上述结果有很大有不同，特别是gamma值= 330.0193）结果如下：
Using maximum autoregressive order for low regime: mL = 3
Using maximum autoregressive order for high regime: mH = 3
Using only first 16 elements of thVar
Performing grid search for starting values...
Starting values fixed: gamma =  100 , th =  0.2919505 ; SSE =  103.1979
Grid search selected lower/upper bound gamma (was:  1 100 ]).
                                          Might try to widen bound with arg: 'starting.control=list(gammaInt=c(1,200))'
Convergence problem code 1. You might want to increase maximum number of iterations by setting 'control=list(maxit=1000)'
Optimized values fixed for regime 2  : gamma =  330.0193 , th =  0.2820525 ; SSE =  90.99445
> summary(ndx.lstar)

Non linear autoregressive model

LSTAR model
Coefficients:
Low regime:
    const1     phi1.1     phi1.2     phi1.3
-0.7398994  1.8526579 -1.2560973  0.7310564
High regime:
    const2     phi2.1     phi2.2     phi2.3
2.7868701 -1.0567428  1.0397095 -0.8612935
Smoothing parameter: gamma = 330
Threshold
Variable: external
Value: 0.2821
Residuals:
     Min       1Q   Median       3Q      Max
-5.95810 -1.28246 -0.17115  1.28210  4.00887
Fit:
residuals variance = 4.789,  AIC = 50, MAPE = 132%
Coefficient(s):
         Estimate  Std. Error  t value  Pr(>|z|)   
const1  -0.739899    0.995576  -0.7432  0.457368   
phi1.1   1.852658    0.257118   7.2055 5.784e-13 ***
phi1.2  -1.256097    0.581746  -2.1592  0.030836 *  
phi1.3   0.731056    0.318841   2.2929  0.021856 *  
const2   2.786870    2.025909   1.3756  0.168941   
phi2.1  -1.056743    0.347791  -3.0384  0.002378 **
phi2.2   1.039710    0.661091   1.5727  0.115784   
phi2.3  -0.861293    0.390391  -2.2062  0.027368 *  
gamma  330.019336  315.129183   1.0473  0.294984   
th       0.282053    0.013939  20.2354 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Non-linearity test of full-order LSTAR model against full-order AR model
F = 0.11978 ; p-value = 0.94616
Threshold
Variable: external>
我的问题是：
    1、这是多变量的LSTAR模型（平滑转换）。为什么在LSTAR中体现不同系列的不稳性问题呢，或者是不需要作平稳性检验？？
    2、我们知道在对作平滑转换模型前，要对基本的回归作非线性检验。这里看不到这一点；是不是不需要这一过程，或含在过程中！
    3、过程中说：Performing grid search for starting values.也即晶格搜索过程；我不知道这一过程的机理，你能不能写出这个过程的程序来；如何选择的！？
    4、这一转换变除X以外的变量z。如是转换变量是Z的滞后项时，又如何办法。是不是程序中也对转换烃量Z及其滞后有选择！如果没有的话，我们要作出是Z的滞后项时，又如何选择呢？
    谢谢！
   
    2、