do stata in mata

mata: mata clear

sysuse auto, clear

reg price weight trunk

mata

b=st_matrix(“e(b)”)’

b

end

******************************

mata: mata clear

sysuse auto, clear

gen constant=1

capt ssc inst tomata

di in red _rc

tomata price constant weight trunk

mata

X=(weight, trunk, constant)

y=price

b=invsym(X’*X)*X’*y

b

//for comparison:

mata stata reg price weight trunk

end

*****************stata mata interact**********************

mata: mata clear

mata:

void underbar(string scalar mat, real scalar thresh)

{

real scalar i

i = rows(st_matrix(mat))

st_matrix(“B”, (st_matrix(mat) :< thresh) :/ st_matrix(mat) :* (J(i,i,

1) – I(i)) )

}

end

sysuse auto,clear

qui reg price trunk weight length turn

mat list e(V)

mata: underbar(“e(V)”, 10)

mat list B

******************************************************

mata: mata clear

sysuse auto, clear

reg price weight trunk

mata

b=st_matrix(“e(b)”)’

b

end

******************************

mata: mata clear

sysuse auto, clear

gen constant=1

capt ssc inst tomata

di in red _rc

tomata price constant weight trunk

mata

X=(weight, trunk, constant)

y=price

b=invsym(X’*X)*X’*y

b

//for comparison:

mata stata reg price weight trunk

end

*****************stata mata interact**********************

mata: mata clear

mata:

void underbar(string scalar mat, real scalar thresh)

{

real scalar i

i = rows(st_matrix(mat))

st_matrix(“B”, (st_matrix(mat) :< thresh) :/ st_matrix(mat) :* (J(i,i,

1) – I(i)) )

}

end

sysuse auto,clear

qui reg price trunk weight length turn

mat list e(V)

mata: underbar(“e(V)”, 10)

mat list B

******************************************************

extract coef,se,t,pv from stata reg

Maarten Buis’s method to extract estimate from stata reg

——-begin correct code —
sysuse auto, clear
reg price mpg rep78 headroom trunk weight
matrix a = vecdiag(e(V))
matrix b = (e(b)\a)’
svmat b
rename b1 beta
gen se = sqrt(b2)
gen t = beta/se
gen p = 2*ttail(63,abs(t))
——— end correct code—–

Or use May’s code

latex hw tempelate

\documentclass[a4paper,12pt,fleqn]{article} %extreport
\usepackage[top=3cm,bottom=1cm,left=2cm,right=2cm]{geometry}
\usepackage{fancyhdr,color,xcolor,dsfont}
\usepackage{listings,ulem,booktabs,graphics,lastpage}
\usepackage{bm,amsmath,amssymb,amsthm,latexsym}
\usepackage{graphicx,fouriernc,hyperref} %fouriernc

\pagestyle{fancy}
\fancyhf{}
\lhead{Macro theory IV:HW1}
\chead{}
\rhead{ d96627003 Chang,JunXiang \quad \thepage/\pageref{LastPage} }

%\usepackage{pifont} fouriernc mathpazo concmath(it can’t use pdf) mathptmx
\def\be{\begin{enumerate}} % Begin Enumerate
\def\ee{\end{enumerate}} % End Enumerate
\def\en{\item} % ENtry (item)
\def\bi{\begin{itemize}} % Begin Itemize
\def\ei{\end{itemize}} % End Itemize

\def\[#1\]{\begin{align*}#1\end{align*}}
\DeclareSymbolFont{letters}{OML}{ccm}{m}{it}
\DeclareMathAlphabet{\mathcal}{OMS}{cmsy}{m}{n}

\renewcommand{\labelitemi}{$\blacktriangleright$}
\newcommand{\nf}{\textrm}
\newcommand{\dsp}{\displaystyle}
\newcommand{\intff}{\int^{\infty}_{0}}
\newcommand{\ert}{e^{-R(t)}}
\newcommand{\ev}{\mathcal{V}}
\newcommand{\dpa}{\partial}
\newcommand{\lw}{\leftarrow}
\newcommand{\llw}{\Leftarrow}
\newcommand{\rw}{\rightarrow}
\newcommand{\rrw}{\Rightarrow}
\newcommand{\ham}{\mathcal{H}}
\newcommand{\lam}{\pounds}
\newcommand{\grl}{\lambda}
\newcommand{\gra}{\alpha}
\newcommand{\grb}{\beta}
\newcommand{\grd}{\delta}
\newcommand{\grt}{\theta}
\newcommand{\limu}{\lim_{t\rw \infty}}
\newcommand{\cov}{\mathsf{Cov}}
\newcommand{\ex}{\mathds{E}}
\newcommand{\bx}{\pmb{X}}
\newcommand{\bgrb}{\pmb{\beta}}
\newcommand{\bgrt}{\pmb{\theta}}
\newcommand{\bA}{\pmb{A}}
\newcommand{\bu}{\pmb{u}}
\newcommand{\bve}{\pmb{\varepsilon}}
\newcommand{\lp}{\mathds{L}}
\newcommand{\I}{\nf{I}}
\newcommand{\iid}{\stackrel{\mathrm{iid}}{\sim}}
\linespread{1.3}
\pagecolor{white}
\color{black}

%=========================here start=======================%
\begin{document}

time format in r axis

x=c(paste(“81m”, 1:12, sep=”"),paste(“82m”, 1:12, sep=”"))
ind=1:24
y=cumsum(rnorm(24,0))
par(lwd=2)
plot(y~ind,type=’b',xaxt=’n',col=4,ylab=”,xlab=”);grid()

axis(1,seq(1,24,2),labels=F,tcl=-0.3)
#down axis
axis(1,seq(2,24,2),labels=F,tcl=-0.5,lwd.ticks=2)
mtext(x[seq(1,24,2)],1,at=seq(1,24,2),padj=3)
#down axis names
mtext(x[seq(2,24,2)],1,at=seq(2,24,2),padj=1)

read multiples files

author:yhliu
source : ptt/statistics

libname ok  ‘c:\SASex’;
Let path=c:\data;

%MACRO aa;
%DO i=1 %TO 13;
DATA a00&i;
INFILE ” &path\a&i..CSV” dlm=’,';
INPUT ID X1-X9;
%END;
%MEND;
%aa;
*Macro PADZERO 在數字前補0, 使其長度固定.
%Macro PADZERO(number,digit);
%LET _l=%LENGTH(&&&number);
%IF &_l<&digit %THEN
%DO _j=&_l+1 %TO &digit;
%LET &number=0&&&number;
%END;
%MEND;
*以下取 alias 的 Macro 做修改.
%Macro RFILES(prefix,begno,endno,digits=0);
%IF &digits<1 %THEN
%LET digits=%LENGTH(&endno);
%DO i=&begno %TO &endno;
%PADZERO(&i,&digits);
%LET _fname=&prefix&i;
DATA &_fname;
INFILE “&Path&_fname..&Ext”;
*      Path 是檔案路徑, Ext 為延伸檔名, 呼叫此 macro
*      前設定之. 或如 alias 原 macro, 用固定字串亦可.
/*
(put your statements for data processing)
*/
RUN;
%END;
%MEND;

%LET Path=C:\;
%LET Ext=dat;
%RFILES(A,1,10,digits=3)

Speed test

%Matlab

m=ones(200000,10);
ret=0;
tic
for i=1:20,
ret=ret+sum(sum(m,1),2);
end
toc

/*Gauss*/
m=ones(200000,10);
ret=0;

starttime=date;
for i (1,20,1);
ret=ret+sumc(sumr(m));
endfor;
endtime=date;
print etstr( ethsec(starttime,endtime) );

/*OX*/
#include
main(){
decl tim,i,m,ret,k;
m = ones(200000,10);
ret =0;

tim = timer();
for (i=0;i<20;++i){
ret=ret+sumc(sumr(m));
}
print(“time lapsed: “, timespan(tim, timer()), “\n”);
}

#R
# File name: loop_rowSums.r

my.loop <- 20
m.dim <- list(nrow = 200000, ncol = 10)
m <- matrix(1, nrow = m.dim$nrow, ncol = m.dim$ncol)
ret <- 0

start <- Sys.time()
for(k in 1 : my.loop){
ret <- ret + sum(rowSums(m))
}
Sys.time() – start
==========================

Speed:

Matlab(0.119)
OX(0.2)
Gauss(0.14)
R(1.625)

http://www.math.ncu.edu.tw/~chenwc/R_note/index.php?item=loop

mySQL

/* IMPORT DATA CSV*/
PROC IMPORT DATAFILE=”C:\DATA\STATADASE17\HOSP17.CSV”
OUT=MYDATA DBMS=CSV REPLACE;
GETNAMES=YES;
RUN;

/*OUTPUT DATA CSV*/
PROC EXPORT DATA=NEW
OUTFILE=’c:\Data\NEW.csv’;

/*CREATE NEW TABLE ORDER BY*/
PROC SQL;
CREATE TABLE NEW AS
SELECT* FROM MYDATA,IDX
WHERE PRSN INV
ORDER BY TAG;
QUIT;

/*USE CASE WHEN TO
GREATE NEW VARIABLES*/
PROC SQL;
CREATE TABLE NEW AS
SELECT * ,
CASE
WHEN PRSN LT 10 THEN ‘ZOO’
WHEN PRSN GE 10 THEN ‘BZOO’
END AS PD,
CASE
WHEN PRSN BETWEEN 1 AND 20 THEN ‘GOOD’
WHEN PRSN BETWEEN 21 AND 30 THEN ‘GOODT’
ELSE ‘NOT_BAD’
END AS DPI
FROM MYDATA;
QUIT;

/*USE IN / NOT IN*/
PROC SQL;
CREATE TABLE NEW AS
SELECT * FROM MYDATA
WHERE PRSN NOT IN
(8,11,14,18,24,27,28,44,47,50,
52,56,57,61,62,64,65,78,80,91,
100,109,110,114,117,
124,126,132,
133,136,138,
143,145,148,
150,155,157,
160,164,166,167,
170,173,175,176,178,179,
180,181,183,184,185,186)
;
QUIT;
/*IN STATA*/
drop if inlist(prsn,8,11,186)

/*R*/
indx=c(8,11,14,18,24,27,28,44,47,50,
52,56,57,61,62,64,65,78,80,91,
100,109,110,114,117,
124,126,132,
133,136,138,
143,145,148,
150,155,157,
160,164,166,167,
170,173,175,176,178,179,
180,181,183,184,185,186)
two=one[!one$PRSN%in%indx,]

/* USE EXCEPT IF TWO TABLES HAVE SAME COLS*/
DATA INDX;
INPUT PRSN @@;
CARDS;
2 3 4 5
7 8 10 12
14 15 19
31 84 92
;
QUIT;

PROC SQL;
CREATE TABLE NEW AS
SELECT PRSN FROM MYDATA
EXCEPT
SELECT PRSN FROM INDX;
QUIT;

PROC PRINT DATA=NEW(OBS=20);
VAR PRSN;
RUN;

R barplot

 x1 <- c(23.2,34.5,76.3,65.8,12.6)
  x2 <- c(15.6,12.4,21.8,20,5.2)
  A <- gl(5,1,5,labels=c("a1","a2","a3","a4","a5"))
  data <- cbind(x1,x2)
  rownames(data) <- levels(A)
  barplot(x1,names.arg=levels(A))
  barplot(t(data),beside=T,ylim=c(0,100),legend.text=colnames(data),
    col=c("grey50","grey80"),ylab="Fréquence")

Simulate price process

http://www.kevinsheppard.com/matlab/Assignment2.html
use R and OX

data management in household data

#stata
cd “c:\data\”
reserve
keep id itm* amt*
outsheet using “c:\data\inc95t.txt”,comma
restore